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Introduction 

High-intensity  focus  ultrasound  (HIFU)  is  gaining  wider  acceptance  in  noninvasive  or 
minimally  invasive  targeting  of  abnormal  tissues  (e.g.  cancer)  for  destruction.  Piezocomposite 
transducer  technology,  especially  for  phased  arrays,  is  providing  high-quality  HIFU  applicators 
with  increased  bandwidth  atfd  reduced  parasitic  cross  coupling  between  the  array  elements.  In 
addition  to  increasing  the  efficacy  of  HIFU  applicators,  these  technological  enhancements 
allow  for  the  use  of  HIFU  arrays  in  imaging  the  target  region  before,  after,  and  intermittently 
during  lesion  formation.  This  leads  to  a  unique  paradigm  of  image-guided  surgery  with  HIFU 
in  which  the  coordinate  systems  for  both  therapy  and  imaging  are  inherently  registered.  This 
project  investigates  the  feasibility  of  using  piezocomposite  phased  arrays  as  dual-mode 
ultrasound  array  (DMUA)  applicators  for  the  noninvasive  treatment  of  primary  breast  cancers. 
Both  therapeutic  and  imaging  capabilities  of  the  dual-mode  arrays  are  investigated  leading  to  a 
real-time  dual-mode  array  system  to  be  used  in  pursuing  in  vivo  animal  experiments  in  the 
future. 

During  the  last  four  years  (original  funding  period  plus  one-year,  no-cost  extension),  we  have 
demonstrated  the  lesion  formation  and  imaging  capabilities  of  the  DMUA.  This  includes  an 
experimental  verification  of  one  of  the  most  important  advantages  of  DMUA  applicators  for 
noninvasive  surgery.  Specifically,  we  have  demonstrated  the  feasibility  of  image-based 
feedback  for  refocusing  in  the  presence  of  strongly  scattering  critical  targets,  e.g.  major  blood 
vessels  or  ribs  or  bony  structures.  This  work,  which  was  presented  at  the  International 
Symposium  on  Biomedical  Imaging  (ISBT04)  as  an  invited  paper,  demonstrates  the  unique 
advantage  of  DMUA  applicators  over  other  forms  of  guidance.  We  have  also  received  an 
invitation  by  Dr.  Ernest  Fellipa  at  the  Riverside  Research  Institute  to  submit  a  journal  paper 
on  the  topic  in  a  special  commemorative  issue  for  the  late  Dr.  Fredrick  Lizzi.  This  paper  is  in 
preparation  and  will  be  submitted  in  January  2006. 

In  addition,  we  have  developed  new  nonlinear  imaging  methods  based  on  the  2nd  order 
Voleterra  filter  (SoVF)  model  that  led  to  significant  enhancement  in  the  image  quality  of 
DMUA  systems.  We  have  also,  investigated  the  detection  and  mapping  of  HlFU-induced 
lesions  using  commercial  diagnostic  scanners. 

Body 

This  report  is  structured  in  accordance  with  the  approved  statement  of  work  (SoW).  In  what 
follows,  we  give  the  tasks  and  subtasks  of  the  approved  SoW  with  each  subtask  followed 
directly  by  what  has  been  accomplished  with  respect  to  it.  For  the  subtasks  planned  for  years  2 
and  3  of  the  grant  period,  they  are  given  here  for  completeness,  but  no  comments  follow  these 
subtasks. 

Task  1.  Thresholds  for  Thermal  Ablation  of  Breast  and  Fatty  Tissue  (Months  1-12): 

a)  Investigation  of  the  intensity/exposure  threshold  curve;  lesion  size  and  characterization  of 
damage  (1  -  3):  Completed  (Year  1  Report) 

b)  Imaging  of  discrete  thermal  lesions  with  therapeutic  arrays  (I  -  3):  Completed  (Year  1 
Report). 

c)  Long-duration  volumetric  ablation  of  porcine  fatty  tissue  (3  -  9):  Completed  (Year  1  Report). 

d)  Imaging  of  volumetric  ablations  using  the  therapeutic  array  and  diagnostic  scanners  (3  -  12): 
Completed  (Year  2  Report). 
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Task  2.  Treatment  Planning  and  Optimization  of  Volumetric  Ablation  with  Phased  Arrays  in  Fatty 
Tissue  (Months  6  -  24): 

a)  Thermal  modeling  based  on  bioheat  equation  and  Saparato  and  Dewey  thermal  dose  integral 
for  damage:  Discrete  lesions  (6-  12):  Completed  (Year  1  Report) 

b)  Thermal  modeling  for  multiple  lesions  with  variable  levels  of  proximity  and  cooling  time 
between  shots  (12-  18):  Completed  (Year  2  Report). 

c)  Optimization  of  multiple-focus  phased  array  patterns  for  simultaneous  placement  multiple 
discrete  lesions  (12  -  24):  Completed.  Simulation  study  results  are  consistent  with  our 
previous  publications  and  expectations  stated  in  Year  2  Report.  The  current  DMUA  prototype 
is  capable  of  producing  up  to  4  high-quality  foci  simultaneously.  This  should  help  reduce  the 
treatment  time  for  volumetric  lesions  by  a  factor  of  4  for  a  conventional  raster  of  focal  points  is 
used. 

Task  3.  Detection  and  Localization  of  Cavitation  Activity  During  Thermal  Lesion  Formation  in 
Breast  and  Fatty  Tissue  (Months  1  -  24): 

a)  Detection  of  subharmonic  activity  in  single-channel  and  beamformed  data  (1  -  12):  Completed 
(Year  1  Report). 

b)  Localization  of  cavitation  from  beamforming  of  multiple  receiving  channels  (12  -  24):  We 
have  the  electronics  and  the  systems  ready  to  perform  this  task.  However,  due  to  the 
malfunction  of  significant  fraction  of  array  elements,  we  are  hesitant  to  drive  the  DMUA 
prototype  at  high  level  to  produce  cavitation  in  tissue.  Since  we  have  not  yet  reached  an 
agreement  with  the  manufacturer,  we  did  not  receive  a  new  prototype  and  this  goal  could  not 
be  pursued.  We  are  confident  that  we  will  be  able  to  test  this  goal  once  the  new  DMUA  is 
received,  but  this  will  have  to  be  done  under  a  future  research  project.  We  note  that  our  earlier 
results  on  overexposure  conditions  in  HIFU  lesion  formation  have  shown  the  feasibility  of 
detection  and  localization  of  cavitation.  We  only  need  to  show  if  the  DMUA  can  localize  the 
onset  of  cavitation  under  normal  exposure  condition. 

c)  Localization  using  time-frequency  and  related  methods  (6  -  18):  Completed  (Year  2  Report) 
Task  4.  Image  Characterization  of  Thermal  Lesions  in  Breast  and  Fatty  Tissue  (Months  1  -  24) 

a)  Characterization  of  grayscale  images  for  discrete  thermal  lesions  (1  -  6):  Completed  (Year  1 
Report). 

b)  Characterization  of  tissue  dependent  parameters  (1  -  12):  Completed  (Year  1  Report: 
(Expected)  negative  results). 

c)  Characterization  of  second  harmonic  imaging  (12  -  24):  Completed  (Year  2  Report). 

d)  Correlation  of  imaging  parameters  with  histologic  characterization  of  tissue  ablations  (12  - 
24):  We  have  obtained  limited  in  vitro  and  in  vivo  lesion  formation  data  using  a  4-MHz  single¬ 
element  transducer  (Courtesy  of  Focus  Surgery,  Inc.,  Indianapolis,  IN).  We  plan  to  use  this  as 
preliminary  data  for  an  NIH  proposal  on  the  subject.  We  are  currently  developing  new 
quantitative  measurements  for  lesion  detection  and  characterization. 

Task  5.  Image-based  Adaptive  Refocusing  of  Therapeutic  Arrays  in  Inhomogeneous  Fatty  Tissue 
(Months  12-24) 

a)  Focusing  with  reference  to  natural  specular  targets  (12  -  18):  Completed  (Year  2  Report). 

b)  Temperature  feedback  (12  -  24):  Completed  (See  Task  6  subtask  d). 
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c)  Discrete  thermal  lesions  as  beacons  (12  -  24):  RF  data  from  single-shot  lesions  in  porcine 
liver  were  analyzed  and  all  normal  lesion  data  sets  clearly  showed  a  strong  localized  echo  from 
the  near  boundary  of  the  lesion  that  can  be  reliably  used  for  focusing.  We  believe  that  we  can 
further  refine  this  idea  by  creating  extremely  small  lesions  amounting  to  point  reflectors. 
Consequently,  the  quality  of  focusing  can  be  significantly  improved. 


Task  6.  Real-time  Dual-mode  Phased  Array  System  for  Volumetric  Thermal  Ablation  of  Breast  and 
Fatty  Tissue  (Months  1  -  36) 

a)  Design  and  fabrication  of  64-channel  receive  system  (1  -  12):  A  64-channel  diplexer  circuit 
allowing  the  use  of  each  array  element  in  transmit  and  receive  has  been  fabricated  and 
tested.  In  addition,  8-to-l  multiplexer  boards  were  designed  and  fabricated  to  allow  data 
collection  from  8  receiving  elements  on  one  analog-to-digital  converter  (A/D).  This  will 
reduce  the  hardware  cost  of  the  system  without  severely  compromising  the  real-time  nature 
of  image  data  collection.  Please  see  d)  below  for  an  example  of  how  this  new  receiver 
system  was  used. 

b)  Design  and  fabrication  of  transmit/receive  control  circuitry  (1  -  18):  Completed  (Year  1 
Report). 

c)  DSP-based  real-time  beamformer  (12  -  30):  During  the  one-year  cost  extension,  we  have 
begun  to  implement  the  real-time  phased  array  controller  on  a  Spartan  III  FPGA  (field 
programmable  gate  array)  from  Xilinx,  Palo  Alto,  CA.  The  new  driver  offers  phase 
resolution  up  to  1.8  0  phase  and  0.01  amplitude  resolution  (compared  to  11.5°  phase  and 
0.07  amplitude  resolution  with  the  current  system). 

d)  Experimental  testing,  calibration,  and  characterization  of  imaging  system  (24  -  36):  Real¬ 
time  data  collection  from  8  channels  to  a  single  A/D  has  been  demonstrated  during  a  9- 
second  heating  experiment.  Not  only  were  we  able  to  demonstrate  beamforming  but  we 
were  also  able  to  extract  temperature  rise  from  the  focus  location.  The  figure  below  shows 
the  estimated  temperature  change  near  the  geometric  focus  of  the  dual-mode  array  (at  100 
mm)  due  to  the  9-second  heating  pulse  from  the  same  array.  The  open  circles  indicate  the 
instants  when  the  data  was  collected  for  temperature  calculations.  As  can  be  seen  from  the 
figure,  noninvasive  temperature  estimates  were  made  at  3  and  6  seconds  during  heating. 

The  heating  pulse  was  interrupted  for  100  ms  for  data  collection  for  each  measurement 
point.  This  includes  the  communications  time  for  the  controlling  computer.  If  a  real-time 
operating  system  were  used,  it  would  have  taken  less  than  2  milliseconds  to  collect  the 
imaging  data  for  temperature  calculation.  The  acquisition  time  can  be  reduced  further  if  we 
dedicate  more  A/Ds  on  the  receiver  side,  i.e.  one  per  4  receiving  elements  instead  of  1  per 
8. 


Key  Research  Accomplishments 

•  First  demonstration  of  a  DMUA  in  forming  lesions  ex  vivo  and  imaging  the  changes  in  tissue 
echogenicity  from  the  lesion  location  and  surrounding  tissues  (Appendix  I). 

•  Demonstrated  the  nonlinear  nature  of  echoes  from  lesion  location  and  the  advantage  of  harmonic 
and  nonlinear  imaging  (based  on  SoVF)  in  improving  the  image  contrast  of  HIFU-induced 
lesions  (Appendix  II). 
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•  First  one-to-one  comparison  of  imaging  performance  of  a  dual-mode  array  with  a  diagnostic 
transducer  (Appendix  III). 

•  First  demonstration  of  image-based  feedback  for  refocusing  the  array  in  the  presence  of  strongly 
scattering  target  (Appendix  IV). 

Reportable  Outcomes 

•  E.S.  Ebbini,  J.  Coad,  and  J.  Bischof,  “Lesion  formation  and  visualization  using  dual-mode 
ultrasound  phased  arrays,”  Proceedings  of  the  2001  IEEE  Ultrasonics  Symposium,  Volume:  2, 

Oct  8-11, 2002,  Page(s):  1351  -  1354. 

•  H.  Yao,  P.  Phupattaranont,  and  E.S.  Ebbini,  “Enhanced  lesion  visualization  in  image-guided 
noninvasive  surgery  with  ultrasound  phased  arrays,”  Proc.  Of  23rd  IEEE  Int.  Conf.  on  Eng.  In 
Medicine  and  Biology,  pp.  2492  -  2495, 2001. 

•  C.  Steidl,  H.  Yao,  P.  Phukpattaranont,  and  E.S.  Ebbini,  “Dual-mode  ultrasound  phased  arrays 
for  noninvasive  surgery:  post-beamforming  image  compounding  algorithms  for  enhanced 
visualization  of  thermal  lesions,”  Biomedical  Imaging,  2002.  Proceedings.  2002  IEEE 
International  Symposium  on  ,  7-10  July  2002,  Page(s):  429  -432. 

•  H.  Yao,  P.  Phukpattaranont,  and  E.S.  Ebbini,  “Nonlinear  methods  for  visualization  of  HIFU- 
induced  lesions,”  Proceedings  of  the  2nd  International  Symposium  on  Therapeutic  Ultrasound,  29 
July  -  1  August,  Seattle,  WA,  Editors:  Andrew,  Crum,  and  Vaezy,  2002,  Page(s):  282  -  289. 

•  H.  Yao,  P.  Phukpattaranont,  and  E.S.  Ebbini,  “Detection  and  mapping  of  thermal  lesions  using 
dual-mode  ultrasound  phased  arrays,”  Proceedings  of  the  2002  IEEE  Ultrasonics  Symposium, 
Volume:  2,  Oct  8-11,  2002,  Page(s):  1435  -1438. 

•  H.  Yao,  P.  Phukpattaranont,  and  E.S.  Ebbini,  “Nonlinear  imaging  methods  for  characterization 
of  HIFU-induced  lesions,”  Proceedings  of  the  SPIE  vol.  4954,  Thermal  Treatment  of  Tissue: 
Energy  Delivery  and  Assessment  II,  Editor:  T.P.  Ryan,  2003,  Page(s):  183-191. 

•  P.  Phukpattaranont  and  Ebbini,  E.S.,  “Post-beamforming  second-order  Volterra  filter  for  pulse-echo 
ultrasonic  imaging,”  IEEE  Trans,  on  Ultrasonics,  Ferroelectrics,  and  Frequency  Control,  vol.  50(8),  pp. 
987-1001,  Aug.  2003. 

•  H.  Yao  and  E.S.  Ebbini,  “Imaging  with  Large-Aperture  Arrays  with  Heterogeneous  Directive 
Elements,”  Proc.  2003  IEEE  Ultrasonics  Symposium,  vol.,  pp.  1243  -  1246. 

•  H.  Yao  and  E.S.  Ebbini,  “Real-time  Monitoring  of  the  Transients  of  HIFU-induced  Lesions,”  Proc. 

2003  IEEE  Ultrasonics  Symposium,  vol.,  pp.  1006  -  1009. 

•  H.  Yao  and  E.S.  Ebbini,  “Dual-mode  Ultrasound  Phased  Arrays  for  Imaging  and  Therapy,” 
Special  Session  on  Advanced  Methods  in  Ultrasound  Imaging,  ISBI  2004,  Arlington,  VA,  April 
2004,  vol.  1,  pp.  24  -  27,  2004.  {Invited  Paper). 

•  Invited  paper  at  the  147th  Meeting  of  the  Acoustical  Society  of  America  in  New  York,  NY,  May  24 
-  28,  2004  (abstract  entitled  Monitoring  of  HIFU-induced  Lesions  Using  Active  Imaging  Methods 
by  Hui  and  Ebbini). 

•  Yao  and  Ebbini,  “Refocusing  of  Dual-mode  arrays  in  the  presence  of  strongly  scattering  objects,” 
Proc.  2004  IEEE  Ultrasonics  Symposium,  vol.  1,  pp.  239  -  242,  2004. 

•  Yao  and  Ebbini,  “Noninvasive  Localized  Ultrasonic  Measurement  of  Tissue  Properties,”  Proc. 

2004  IEEE  Ultrasonics  Symposium,  vol.  1,  pp.  724  -  727,  2004. 
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•  H.  Yao,  Active  Ultrasonic  Imaging  Methods  for  HIFU-induced  Lesions,  Ph.D.  Dissertation, 
University  of  Minnesota,  2006  (defended  2005). 


Conclusions 

After  characterizing  and  optimizing  the  imaging  capabilities  of  the  DMUA  prototype  in  year  2,  we 
have  experimentally  demonstrated  one  of  the  most  important  advantages  of  these  arrays  in  noninvasive 
surgery.  Namely,  we  have  demonstrated  the  feasibility  of  using  image-based  feedback  to  refocus 
DMUAs  in  the  presence  of  strongly  scattering  targets.  We  have  also  begun  the  investigation  of 
quantitative  imaging  of  HIFU-induced  lesions  through  active  imaging  methods.  Specifically,  we  are 
developing  techniques  for  measuring  local  absorption,  perfusion,  and  elastic  properties  of  HIFU- 
induced  lesions.  We  are  developing  these  techniques  using  both  diagnostic  scanners  and  our  DMUA 
prototypes. 
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Abstract-  A  1  MHz  spherical-section  64-element  linear  plezo- 
composlte  phased  array  and  a  supporting  driving  and  data  ac¬ 
quisition  system  were  recently  tested  for  use  as  a  dual-mode 
HIFU  applicator  system.  In  therapeutic  mode,  the  array  was 
shown  to  be  capable  of  producing  focal  Intensities  In  excess  of 
3500  W/cm2  at  its  geometric  center  (100  mm  radius  of  cur¬ 
vature).  Imaging  tests  of  quality  assurance  tissue-mimicking 
phantoms  as  well  as  computer  simulations  confirmed  that  the 
array  has  an  oval-shaped  Imaging  field  of  view  ( IxFOV)  with 
nearly  50  dB  dynamic  range  extending  from  70  to  120  mm 
in  the  axial  direction  and  ±  20  mm  in  the  transverse  direc¬ 
tion  [1].  This  is  larger  than  its  therapeutic  operating  field 
(ThxOF),  which  Is  defined  as  the  region  around  the  geometric 
focus  where  the  loss  In  intensity  gain  is  within  1  dB  from  the 
Intensity  gain  at  the  geometric  center.  We  have  tested  the  ar¬ 
ray  imaging  capability  in  the  visualization  of  the  formation  of 
discrete  and  volumetric  lesions  in  freshly  excised  porcine  liver 
samples.  Experimental  results  clearly  show  that  the  enhance¬ 
ment  of  the  tissue  echogenecity  at  the  lesion  location  Is  more 
pronounced  In  the  harmonic  Images  when  compared  to  images 
at  the  fundamental.  Results  also  clearly  show  that  the  har¬ 
monic  imaging  mode  produces  more  accurate  mapping  of  the 
lesion  size  and  shape  as  determined  by  histologic  evaluation. 

1.  INTRODUCTION 

Image  guidance  has  long  been  recognized  as  the  “enabling 
technology”  for  noninvasive  or  minimally  invasive  thermal 
surgery.  Highly  refined  energy  application  devices  have 
been  introduced  for  many  years.  However,  without  reli¬ 
able  imaging  techniques  for  visualization  of  thermal  lesions, 
noninvasive  thermal  surgery  has  failed  to  find  widespread 
acceptance  in  the  clinic.  Recently,  image  guidance  methods 
based  on  well  established  imaging  modalities  like  MRI[8], 
CT[9],  or  ultrasound  [10,  5]  have  been  proposed.  Other 
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imaging  modalities  are  also  being  developed  and  may  be¬ 
come  available  in  the  foreseeable  future  [2], 

An  area  unique  to  ultrasound  that  could  revolutionize 
the  field  of  image  guided  surgery  is  the  development  of 
a  new  generation  of  dual-mode  high-power  phased  array 
systems  capable  of  both  imaging  and  therapy  [1],  These 
piezocomposite  transducers  can  produce  focal  intensity  lev¬ 
els  needed  for  ablative  and  coagulative  thermal  surgery  with 
high  precision.  Furthermore,  the  operating  bandwidth  of 
such  transducers  allows  for  imaging  the  treatment  region 
with  adequate  image  quality  to  delineate  important  land¬ 
marks  within  and  around  the  target  volume.  With  these  ca¬ 
pabilities,  it  is  possible  to  operate  these  arrays  in  a  “self¬ 
registration”  mode  whereby  the  array’s  imaging  capabilities 
are  utilized  in  characterization  of  the  tissue  response  pre¬ 
cisely  at  the  expected  lesion  location.  This  is  due  to  the  fact 
that  the  beamforming  is  common  to  both  the  imaging  and 
therapy  modes. 

In  this  paper  we  present  results  of  lesion  formation  and  vi¬ 
sualization  using  a  64-element  1  MHz  array  originally  de¬ 
signed  for  therapeutic  applications  [4], 

2.  DUAL-MODE  PHASED  ARRAY  SYSTEM  FOR 
NONINVASIVE  SURGERY 

A  variety  of  nonplanar  phased  array  geometries  for  thera¬ 
peutic  applications  have  been  investigated  thoroughly  since 
the  late  1980s  [7].  Briefly,  optimally  designed  therapeutic 
arrays  have  the  following  geometrical  characteristics: 

•  Nonplanar  with  radius  of  curvature  and  /-number  de¬ 
termined  by  the  depth  and  volume  of  the  ThxOF  (/- 
numbers  typically  ~  1). 

•  Coarsely  sampled  with  element  dimensions  ~  1  -  3A 
to  maintain  high  efficiency.  Maximum  element  size 
is  limited  by  the  extent  of  ThxOF ;  element  directivity 
variation  within  ThxOF  is  ~  1  dB. 
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Unfortunately,  these  characteristics  render  the  imaging  ca¬ 
pabilities  of  therapeutic  arrays  poor  due  to  their  unaccept¬ 
ably  high  grating  lobe  levels.  However,  it  is  possible  to  de¬ 
sign  nonplanar  therapeutic  arrays  with  limited  FOV  cover- 
ing  or  exceeding  the  extent  of  their  ThxOFs.  This  limited 
imaging  capability  can  be  quite  significant  in  that  it  poten¬ 
tially  allows: 

•  Visualizing  critical  structures  in  the  vicinity  of  the  tar¬ 
get  point(s),  e.g.,  blood  vessels,  bone  or  gas  spaces. 

•  Measuring  temperature  dependent  changes  that  may 
allow  the  assessment  of  the  quality  of  the  therapeutic 
beam  in  situ  at  sub-therapeutic  (i.e.,  nondestructive) 
levels. 

•  Monitoring  the  effects  of  the  therapeutic  beam  during 
and/or  immediately  after  lesion  formation. 

We  have  designed  and  tested  a  64-element  ID  array  on  a 
spherical  surface  as  a  prototype  for  a  dual-mode  transducer 
for  noninvasive  surgery.  In  previous  publications  [1,  5,  6], 
we  have  demonstrated  the  imaging  capabilities  of  this  ar¬ 
ray  and  its  IxFOV.  In  addition,  we  have  shown  that  it  is 
possible  to  perform  noninvasive  temperature  imaging  using 
this  array  [1]  .In  this  report,  we  demonstrate  the  feasibility  of 
imaging  thermal  lesions  in  freshly  excised  tissue  samples. 

2.1.  Imaging  Characterization 

2.1.1.  Data  collection  system 

The  data  system  collection  for  purposes  of  image  forma¬ 
tion  is  shown  in  Figure  1.  The  array  is  connected  to  a  64- 
channel  pulser/phased  array  driver  capable  of  driving  the 
array  in  continuous  or  pulsed  mode  with  time  resolution  of 
25  nanoseconds.  Channel  delays  from  25  ns  to  1  ms  can 
be  realized  with  this  driver.  In  its  current  form,  the  receiver 
channel  data  is  connected  to  a  23-bit  20  MSample/second 
digitizer  through  a  4  x  64  matrix  switch. 

2.1.2.  Image  Formation  Model 

Figure  1  summarizes  the  image  acquisition  and  image  for¬ 
mation  model.  A  64-element  array  optimized  for  maximum 
energy  delivery  at  1  MHz  operating  frequency  is  used  for 
lesion  formation  in  sample  tissue.  Lesions  are  formed  by 
focusing  the  array  at  a  point  within  the  target  and  main¬ 
taining  high-power  output  for  time  intervals  on  the  order 
of  seconds  (1-5  seconds  typical).  The  power  is  inter¬ 
rupted  for  short  intervals  (milliseconds)  to  acquire  image 
data  by  transmitting  short  (fts)  pulses  from  all  64-elements 
and  receiving  on  selected  elements  using  a  matrix  switch. 
Once  the  image  data  set  is  collected,  RF  beamforming  is 
performed  to  form  standard  echographic  images  of  the  tar¬ 
get  region.  It  is  well  known  that  these  images  have  very 
low  contrast  due  to  the  coherent  nature  of  image  acquisition 
which  produces  speckle.  This  renders  the  lesion  detection 


Fig.  1.  Experiment  setup  and  image  reconstruction  algorithm  for 
the  64-e!ement  dual-mode  array. 


very  difficult  using  ultrasound.  However,  due  to  the  nonlin¬ 
earity  of  echoes  from  freshly  formed  thermal  lesions,  (sec¬ 
ond)  harmonic  imaging  provides  significant  improvement  in 
contrast  for  purposes  of  lesion  detection  and  mapping. 

2.1.3.  Harmonic  Imaging 

Due  to  the  double  resonance  characteristics  of  the  64-element 
array  described  in  [1],  it  is  easy  to  form  second  harmonic 
images  by  post-beamforming  filtering  of  RF  data  to  elimi¬ 
nate  the  fundamental.  Of  course,  second  harmonic  imaging 
requires  the  full  transmit  power  be  used  for  transmit  beams 
as  opposed  to  the  synthetic  aperture  approach  in  which  each 
array  element  is  driven  separately  as  a  transmitter.  There¬ 
fore,  for  harmonic  imaging,  we  have  used  the  64-channel 
driver  to  form  full  transmit  beams  directly  for  each  recep¬ 
tion.  The  images  obtained  using  this  technique  are  referred 
to  as  single-transmit  images.  These  images  are  equivalent 
to  full  data  sets  along  the  transmit  focus  direction,  but  are 
otherwise  compromised,  While  this  may  seem  like  an  un¬ 
acceptable  degradation,  it  does  have  distinct  advantages  in 
image-guided  surgery.  When  the  transmit  imaging  beam 
has  the  same  structure  as  the  therapeutic  beam,  single-transmit 
images  allow  for  the  visualization  of  any  critical  structures 
near  or  along  the  line  of  sight  of  the  target  point.  This  can 
be  used  to  modify  the  focusing  and  steering  parameters  of 
the  therapeutic  beam  to  avoid  or  target  these  critical  struc¬ 
tures.  It  is  also  useful  in  isolating  specular  reflectors,  e.g., 
large  blood  vessels,  that  may  be  useful  in  characterizing  the 
quality  of  the  therapeutic  beam  in  the  vicinity  of  the  target. 
As  we  show  below,  this  mode  is  extremely  useful  in  visu¬ 
alization  of  discrete  thermal  lesions  when  used  in  both  the 
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fundamental  and  in  the  second  harmonic  modes. 

2.2.  Therapeutic  Characterization 

The  64-element  piezocomposite  array  described  above  was 
extensively  tested  as  a  high-power  therapeutic  applicator. 
Discrete  and  volumetric  thermal  lesions  were  formed  in- 
vitro  bovine  muscle  [4],  More  recently,  we  have  performed 
a  series  of  experiments  for  forming  discrete  and  volumet¬ 
ric  lesions  in  porcine  liver  using  the  setup  shown  in  Fig¬ 
ure  1,  To  illustrate  one  of  the  unique  capabilities  of  phased 
arrays  in  forming  lesions  of  various  shapes,  Figure  2  shows 
two  cigar-shaped  lesions  formed  in  porcine  liver  simulta¬ 
neously  using  a  double-focus  pattern  with  focal  intensity  of 
660  W/cm2  at  each  focus  applied  for  4  seconds.  It  has  been 
determined  that  this  array  is  capable  of  producing  focal  in¬ 
tensity  levels  of  over  3500  W/cm2  in  well-degassed  water. 
This  level  of  focal  intensity  is  sufficient  to  produce  a  full 
range  of  known  therapeutic  effects  of  ultrasound, 


Fig.  2.  Two  cigar-shaped  lesions  formed  in  porcine  liver  sample 
after  4-second  exposure  to  a  double-focus  beam  pattern  produced 
by  (he  64-element  therapeutic  array. 


3.  HISTOLOGIC  EVALUATION 

The  liver  was  cut  with  approximate  square  cross  section  of 
4.5  x  4.5cm2  so  they  fit  snugly  within  the  sample  holder 
without  being  squeezed.  Samples  were  brought  in  contain¬ 
ers  in  saline  buffer  within  an  hour  of  sacrificing  the  animal. 
Histological  measurements  were  made  for  identified  lesion 
and  comparisons  with  gross  measurements  were  made.  Cross 
sectional  cuts  into  the  samples  were  made  in  planes  parallel 
to  the  face  of  the  array  (i.e.,  at  fixed  axial  distance)  to  insure 
finding  each  lesions  that  is  formed.  The  following  zones  of 
damage  were  identified: 

•  At  the  Center  of  each  lesion,  cells  undergo  vital  en¬ 
zyme  protein  denaturation.  For  higher  dose  levels, 
ablation  and  evidence  of  water  vaporization  can  be 
observed  in  this  central  zone  as  well. 

•  An  outer  ring  of  pyknotic  nuclei  with  thickness  be¬ 
tween  0.1  -  0.5  mm  is  variably  present  for  all  lesion. 

•  An  inner  ring  with  nuclei  in  this  zone  show  a  transi¬ 
tion  between  pyknotic  nuclei  of  the  outer  ring  and  the 


central  zone  cells  that  have  undergone  vital  enzyme 
protein  denaturation. 

A  variety  of  exposure  levels  and  durations  varying  from  un¬ 
der  exposure  (low  probability  of  lesion  formation)  to  over 
exposure  (strong  cavitation  and  distorted  lesion  formation) 
were  examined  at  histological  cross  sections  2  mm  apart 
starting  5  mm  into  the  tissue  sample  covering  the  full  ex¬ 
tent  of  the  lesion.  The  far  end  of  the  lesion  was  always  a 
point  and  it  occurred  at  the  geometric  center  of  the  array. 
The  near  end  and  the  maximum  lateral  cross  section  of  the 
lesion  depend  on  the  exposure  level.  Normal  exposure  lev¬ 
els  produced  expected  cigar-shaped  lesions  of  length  pro¬ 
portional  to  intensity.  Over  exposures,  on  the  other  hand, 
produced  inverted-cone  lesions  with  the  near  end  extending 
to  the  front  surface  of  the  sample  in  some  cases  (surface 
bums). 

4.  LESION  VISUALIZATION 

Figure  3  shows  25-dB  gray  scale  images  obtained  using  the 
single-transmit  imaging  mode  (transmit  focus  at  0  lateral 
and  100  axial)  obtained  before  (left)  and  after  (right)  sin¬ 
gle  lesion  formation  (4  seconds  at  850  W /cm2).  The  top 
pair  was  formed  at  the  fundamental  and  the  bottom  pair  was 
formed  at  the  second  harmonic.  For  each  pair,  the  before 
and  after  images  were  normalized  with  respect  to  the  maxi¬ 
mum  in  the  image  after. 

The  images  acquired  after  lesion  formation  both  show  in¬ 
creased  echogenecity  at  the  lesion  location  starting  at  92 
mm.  This  is  consistent  with  the  histologic  evaluation  of  the 
sample.  The  images  also  suggest  that  the  extent  of  the  le¬ 
sion  is  8  mm,  also  in  accordance  with  histologic  evaluation. 
It  should  be  mentioned,  however,  that  the  boundaries  of  the 
lesion  are  much  easier  to  detect  in  the  second  harmonic  im¬ 
age. 

Figure  4  shows  25-dB  gray  scale  images  obtained  using  the 
single-transmit  imaging  mode  (transmit  focus  at  0  lateral 
and  100  axial)  obtained  before  (left)  and  after  (right)  single 
lesion  formation  (4  seconds  at  1200  W /cm2).  This  is  only  a 
20%  change  in  dose  delivered  to  the  focus  compared  to  the 
previous  case.  Here  the  images  acquired  after  lesion  forma¬ 
tion  suggest  a  lesion  of  10  mm  length  starting  at  90  mm, 
both  in  accordance  with  histological  evaluation.  As  in  the 
previous  case,  the  second  harmonic  images  show  significant 
increase  in  echogenecity  at  the  lesion  location  compared  to 
the  fundamental. 

5.  CONCLUSIONS 

The  first  experiment  represent  a  fairly  homogeneous  target 
tissue  while  the  second  has  a  pair  of  blood  vessels  at  the 
target  location  as  can  be  seen  from  the  images  taken  be¬ 
fore  lesion  formation.  While  both  the  fundamental  and  har¬ 
monic  images  show  enhanced  lesion  contrast  in  both  cases, 
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Fig.  3.  Fundamental  (top)  and  second  harmonic  (bottom)  images 
before  (left)  and  after  lesion  formation  using  a  single  focus  (4  sec¬ 
onds  at  850  W/cm2). 


Fig.  4.  Fundamental  (top)  and  second  harmonic  (bottom)  images 
before  (left)  and  after  lesion  formation  using  a  single  focus  (4  sec¬ 
onds  at  1200  W/cm2). 


the  harmonic  images  show  a  net  increase  in  contrast  by  22 
dB  in  the  second  case.  On  the  other  hand,  the  net  increase 
in  contrast  at  the  fundamental  is  about  7  dB.  In  addition,  the 
spatial  definition  of  the  lesion  in  the  harmonic  images  is  su¬ 
perior  to  the  fundamental. 

The  results  shown,  in  addition  to  numerous  similar  experi¬ 
ments,  suggest  that  the  second  harmonic  single-transmit  im¬ 
ages  can  be  used  to  detect  even  small  changes  in  the  start¬ 
ing  point  of  formed  lesions.  This  applies  to  a  variety  of 
exposure  conditions,  including  over  exposure.  The  results 
of  these  experiments  also  show  that  the  lesion  images  are  in 
good  agreement  with  histological  examinations,  including 
lesion  shapes.  These  results  will  be  reported  in  an  upcom¬ 
ing  publication. 
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APPENDIX  II 

Post-Beamforming  Second-Order  Volterra 
Filter  for  Pulse-Echo  Ultrasonic  Imaging 
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Abstract — We  present  a  new  algorithm  for  deriving  a 
second-order  Volterra  filter  (SVF)  capable  of  separating  lin¬ 
ear  and  quadratic  components  from  echo  signals.  Images 
based  on  the  quadratic  components  are  shown  to  provide 
contrast  enhancement  between  tissue  and  ultrasound  con¬ 
trast  agents  (UCAs)  without  loss  in  spatial  resolution.  It 
is  also  shown  that  the  quadratic  images  preserve  the  low 
scattering  regions  due  to  their  high  dynamic  range  when 
compared  with  standard  B-mode  or  harmonic  images.  A 
robust  algorithm  for  deriving  the  filter  has  been  developed 
and  tested  on  real-time  imaging  data  from  contrast  and 
tissue-mimicking  media.  Illustrative  examples  from  image 
targets  containing  contrast  agent  and  tissue-mimicking  me¬ 
dia  are  presented  and  discussed.  Quantitative  assessment 
of  the  contrast  enhancement  is  performed  on  both  the  RF 
data  and  the  envelope-detected  log-compressed  image  data. 
It  is  shown  that  the  quadratic  images  offer  levels  of  enhance¬ 
ment  comparable  or  exceeding  those  from  harmonic  filters 
while  maintaining  the  visibility  of  low  scattering  regions  of 
the  image. 


I.  Introduction 

Increasing  interest  in  extending  the  capabilities  of  ul¬ 
trasound  imaging  by  utilizing  ultrasound  contrast  agents 
(UCAs)  has  heightened  the  need  for  more  suitable  imaging 
techniques.  In  standard  B-mode  imaging,  UCAs  increase 
the  echogenicity  from  perfused  tissues  [1].  This  results  in 
improved  endocardial  border  detection  in  left  ventricular 
opacification,  which  leads  to  a  better  analysis  of  wall  mo¬ 
tion  abnormalities  [2].  Nevertheless,  in  the  myocardium 
where  the  ratio  of  blood  volume  to  tissue  is  quite  low 
(approximately  10%  [3]),  the  backscatter  from  the  small 
number  of  microbubbles  in  vessels  can  be  dominated  by 
echoes  from  surrounding  tissue.  In  this  occurrence,  stan¬ 
dard  B-mode  imaging  offers  inferior  UCA  detectability  in 
the  presence  of  tissue,  stated  as  agent-to-tissue  ratio  [4]. 
In  order  to  increase  the  sensitivity  of  UCA  detections,  var¬ 
ious  new  imaging  techniques  [5],  [6]  have  been  developed 
by  employing  some  specific  acoustic  signatures  of  UCAs, 
such  as  nonlinear  and  transient  scattering. 

Imaging  techniques  based  on  nonlinear  oscillations  have 
been  designed  for  separating  and  enhancing  nonlinear 
UCA  echoes  from  a  specified  region  of  interest  within  the 
imaging  field,  including  second  harmonic  (SH)  B-mode 
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imaging  and  pulse  inversion  (PI)  Doppler  imaging  [7]. 
The  SH  imaging  employs  a  fundamental  frequency  trans¬ 
mit  pulse  and  produces  images  from  the  second  harmonic 
component  of  received  echoes  by  using  a  second  harmonic 
bandpass  filter  (BPF)  to  remove  the  fundamental  fre¬ 
quency.  In  order  to  increase  UCA  detection  sensitivity  in 
the  limited  transducer  bandwidth  condition,  spectral  over¬ 
lap  between  fundamental  and  second  harmonic  parts  need 
to  be  minimized  by  transmitting  narrow-band  pulses  re¬ 
sulting  in  an  inherent  tradeoff  between  contrast  and  spa¬ 
tial  resolution. 

In  PI  imaging,  a  sequence  of  two  inverted  acoustic 
pulses  with  appropriate  delay  is  transmitted  into  tissue. 
Images  are  produced  by  summing  the  corresponding  two 
backscattered  signals.  In  the  absence  of  tissue  motion,  the 
resulting  sum  can  be  shown  to  contain  only  even  harmonics 
of  the  nonlinear  echoes  [7].  The  PI  imaging  overcomes  the 
tradeoff  between  contrast  and  spatial  resolution  because  it 
utilizes  the  entire  bandwidth  of  the  backscattered  signals 
[7].  As  a  result,  superior  spatial  resolution  can  be  achieved 
when  compared  with  SH  imaging.  Moreover,  it  has  been 
shown  that  PI  imaging  can  be  operated  in  a  continuous 
imaging  mode  with  low  mechanical  indices  (Mis)  [8].  The 
PI  imaging  is  sensitive  to  tissue  motion  because  it  is  a 
multiple  pulse  technique;  therefore,  PI  detection  is  com¬ 
bined  with  Doppler  detection  leading  to  a  new  technique 
called  PI  Doppler.  The  PI  Doppler  utilizes  the  advantages 
from  both  detection  schemes  and  circumvents  the  tissue 
motion  problem  [7].  Nevertheless,  an  inherent  multipulse 
technique  of  PI  imaging  results  in  the  reduction  of  imaging 
frame  rates. 

In  order  to  detect  backscattered  signals  due  to  UCAs  in 
pulse-echo  ultrasound  imaging  with  single  transmit  pulse 
per  line  and  overcome  the  tradeoff  between  contrast  and 
spatial  resolution,  we  have  developed  a  new  imaging  tech¬ 
nique  based  on  the  Volterra  filter  [9].  The  Volterra  filter 
is  a  dynamic  filter  that  operates  in  parallel  on  the  lin¬ 
ear,  quadratic,  cubic,  etc.,  signal  components  to  produce 
its  output.  It  has  a  discrete  convolutional  form  with  finite 
memory  that,  is  commonly  used  in  nonlinear  digital  signal 
processing  (DSP)  applications  [10].  While  the  Volterra  fil¬ 
ter  output  is  nonlinear  in  terms  of  the  input  data,  it  is 
linear  in  terms  of  its  coefficients.  As  a  result,  linear  pro¬ 
cessing  can  be  applied  to  identify  Volterra  kernels  (i.e.,  fil¬ 
ter  coefficients).  The  identification  of  Volterra  kernels  has 
been  demonstrated  in  several  applications  [10],  [11].  For 
example,  the  linear  and  quadratic  Volterra  kernels  identi¬ 
fied  by  an  adaptive  filtering  algorithm  based  on  recursive 
least-squares  approach  of  a  second-order  Volterra  model 


08S5-3010/S10.00  ©  2003  IEEE 


988 


IEEE  TRANSACTIONS  ON  ULTRASONICS,  FERROELECTRICS,  AND  FREQUENCY  CONTROL,  VOL.  50,  NO.  8,  AUGUST  2003 


are  applied  to  separate  the  linear  and  quadratic  responses 
of  a  tension  leg  platform  [12]. 

In  this  paper,  we  present  postbeamforming  nonlinear 
filter  based  on  the  second-order  Volterra  filter  (SVF).  The 
filter  is  capable  of  separating  linear  and  quadratic  compo¬ 
nents  from  UCA-backscattered  signals.  Filter  coefficients 
c.a.n  be  identified  by  forming  a  system  of  linear  equations 
from  a  beamformed  RF  data  segment  in  the  UCA  or  tis¬ 
sue  region  and  obtained  by  solving  a  minimum-norm  least- 
squares  (MNLS)  problem.  It  is  shown  that  the  system  iden¬ 
tification  approach  to  a  nonlinear  ultrasonic  imaging  is  ro¬ 
bust  and  can  be  applied  to  form  nonlinear  images  through¬ 
out  the  imaging  field,  i.e.,  not  confined  to  the  region  where 
the  filter  coefficients  were  derived.  Images  produced  from 
quadratic  output  of  the  SVF  model  exhibit  excellent  con¬ 
trast  enhancement  without  loss  in  spatial  resolution.  In 
addition,  quadratic,  images  have  increased  dynamic  range 
compared  to  standard  B-mode  and  harmonic  images  allow¬ 
ing  the  preservation  of  image  features  along  with  contrast 
improvement. 


II.  Theory 

For  simplicity  and  without  loss  of  generality,  we  illus¬ 
trate  the  nonlinear  postbeamforming  algorithm  using  the 
SVF  model.  The  algorithm  described  in  this  section  for  de¬ 
riving  the  filter  coefficients  of  the  SVF  extends  to  higher 
order  in  a  straightforward  manner.  Initial  experience  with 
this  model  indicate  that  it  is  largely  sufficient  for  tissue  re¬ 
sponse.  Furthermore,  the  SVF  is  computationally  feasible 
for  real-time  application  with  today’s  technology,  which 
makes  it  more  attractive  from  the  implementation  point 
of  view. 

A.  Second- Order  Volterra  Model 

Results  from  [9]  have  shown  the  validity  of  a  SVF  as 
a  model  for  pulse-echo  ultrasound  pulse-echo  data  from 
tissue  mimicking  media.  An  input-output,  system  identifi¬ 
cation  approach  was  used  to  estimate  the  coefficients  of 
the  linear  and  quadratic  components  of  the  SVF  in  the 
frequency  domain.  However,  while  the  system  identifica¬ 
tion  study  was  necessary  to  establish  the  applicability  of 
SVF  to  ultrasound  pulse-echo  data,  it  is  not  useful  for 
imaging  purposes  as  it  requires  access  to  both  the  input 
and  the  echo  data  from  distinct  scatterers  in  the  tissue- 
mimicking  media.  An  appropriate  approach  for  imaging 
operates  on  the  beamformed  RF  data  to  separate  the  lin¬ 
ear  and  quadratic  components  regardless  of  the  input.  This 
signal  separation  approach  allows  us  to  extract  the  linear 
and  quadratic  signal  components  from  the  beamformed 
data  to  form  linear  and/or  quadratic  images  separately  or 
compounded.  In  this  paper,  we  emphasize  the  quadratic 
images  obtained  from  the  SVF  model  based  on  a  linear 
and  quadratic  prediction  model  of  beamformer  output  de¬ 
scribed  in  the  remainder  of  this  section. 


uL(n) 


Fig.  1.  Separation  of  beamformed  R.F  data  into  linear  and  quadratic 
components  using  the  SVF. 


1.  Signal  Separation  Model:  Fig.  1  shows  a  simple  block 
diagram  of  the  imaging  system  based  on  SVF.  The  SVF 
operates  on  the  beamformer  output  to  produce  the  linear 
and  quadratic  components,  and  uq(ii),  respectively. 
Estimates  of  the  total  beamformer  output  can  be  obtained 
from  these  components  simply  by  adding  them 

u{n)  =  Mn)  +  uq(u),  (1) 

where  u(n),  Ui(n),  and  UQ{n)  are  the  total,  linear,  and 
quadratic  estimations,  respectively.  The  separation  of  the 
linear  and  quadratic  components  can  be  achieved  once 
the  coefficients  of  the  kernels,  /i£,(i)  (linear),  and  hQ(j,k) 
(quadratic)  are  found.  In  the  following  subsection,  we  de¬ 
scribe  a  MNLS  approach  for  determining  these  coefficients. 

2.  The  MNLS  Estimation  of  SVF  Coefficients:  The  re¬ 
sponse  of  a.  quadratically  nonlinear  system  with  memory, 
u(n  +  1)  can  be  predicted  by  a  (discrete)  second-order 
Volterra  model  operating  on  the  m  past  samples  as  fol¬ 
lows 


m  —  1 

u(n  4- 1)  =  ^  u{n  —  i)hi(i) 
i=0 

m— 1 771 —  1 

+  Y!  u(n  "  J'Mn  “  k)hQti> k)>  (2) 

j= o  k=j 

where  hffii)  is  linear  filter  coefficients,  and  fiQ(j,k)  repre¬ 
sents  quadratic  filter  coefficients.  Note  that  while  u(n+  1) 
is  nonlinear  with  respect  to  the  beamformed  data,  it  is 
linear  with  respect  to  the  coefficients  of  the  linear  and 
quadratic  kernels  of  the  SVF.  Recognizing  this  fact,  one 
can  rewrite  (2)  in  vector  form 

u(n  +  1)  =  uT(n)h,  (3) 

where  the  data  vector,  u(n),  is  defined  at  sample  n  as 

u(n)  =  [tt(n),  u{n  —  1),  u(n  —  2), . . .  ,  u(n  —  m  +  1), 
u2(n),  u{n)u{n  —  1), ...  ,  u2(n  —  m  +  l)]r, 

and  the  filter  coefficient  vector,  h,  can  be  expressed  as 

h  =  [M°)>  Ml),  M2),  ■  ■  ■  ,  hL(m  -  1), 

M°<°),  hQ(0,l),...  ,  hQ (m  —  l,m  —  1)]T, 
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where  m  is  the  system  order  and  superscript  T  denotes  the 
transpose.  The  total  number  of  independent  filter  coeffi¬ 
cients,  N,  is  equal  to  (m2  +  3m)/2  assuming  a  symmetri¬ 
cal  quadratic  kernel  (i.e.,  h.Q(j,k)  =  h,Q(k,j)).  Similarly, 
u{n  +  2),  u{n  +  3), . . . ,  u{n  +  M)  can  be  represented  in  the 
form  of  (3)  and  expressed  in  the  matrix  form 

f  =  Gh,  (4) 

where  the  vector  f  and  the  matrix  G  are  defined  as 
f  =  [«(n  +  1),  u(n  +  2),  . . . ,  u(n  +  M)]r 
and 


B.  Regularization 

The  SVD  of  G  forms  a  basis  for  regularization  by  ap¬ 
propriate  selection  of  singular  modes  that  enhance  the  re¬ 
constructed  image  in  some  sense.  There  are  a  number  of 
approaches  for  regularization  of  (8),  including  single  pa¬ 
rameter  and  rank  reduction  regularization  [14].  The  lat¬ 
ter,  sometimes  referred  to  as  the  truncated  singular  value 
decomposition  (TSVD),  produces  a  solution  by  truncat¬ 
ing  the  number  of  singular  modes  of  G  with  the  smallest 
singular  values  below  a  certain  threshold.  The  fcth  order 
TSVD  solution  is  given  by 


G  =  [u(n),  u (n  +  1), . . . ,  u (n  +  M  -  1)]T, 

where  M  is  the  number  of  linear  equations  (observations). 
The  linear  and  quadratic  filter  coefficients  can  be  esti¬ 
mated  by  seeking  an  appropriate  solution  of  (4).  Well- 
known  solutions  to  (4)  are  the  least  squares  (LS)  solu¬ 
tion  for  the  overdetermined  case  (more  constraints  than 
unknowns)  and  the  minimum  norm  (MN)  solution  in  the 
underdetermined  case  (less  constraints  than  unknowns).  A 
MNLS  solution  can  be  obtained  by 

hMNLS  =  G*f,  (5) 


h*  =  £ 


(9) 


where  the  truncation  parameter  k  <  r,  also  known  as  the 
rank  of  the  approximation,  is  the  number  of  singular  modes 
used  to  compute  the  estimate. 

The  regularization  is  guided  by  the  mean  square  error 
criterion,  approximated  by 

m  =  loiogio  ('-f--^hfc  l|2) ,  (io) 


where  G  f  is  a  generalized  inverse  [13].  This  solution  ap¬ 
plies  to  both  the  overdetermined  and  underdetermined 
cases  as  it  accounts  for  the  effective  rank  of  the  matrix 
G.  In  general,  G  could  be  rank  deficient,  i.e.,  has  rank 
v  <  min {M,N}.  When  this  is  the  case,  there  are  infinite 
number  of  solutions  to  (4)  that  produce  the  same  LS  er¬ 
ror,  ||f  —  GhLs||2.  The  MNLS  is  the  unique  solution  to 
(4)  with  minimum  norm,  i.e.,  ||hMNLs||2  <  ||h£,s||2  Vh ls- 
The  singular  value  decomposition  (SVD)  of  G  is  given  by 


G  =  U£VT 

=  ^2<?i  u<vf, 

i=l 


(6) 


and 


Gf  =  VEfUr 


(7) 


where  £  is  a  M  x  N  diagonal  matrix  with  singular  val¬ 
ues  o\  >  02  >  ■  ■  ■  >  crr  >  ov+i  =  . . .  =  (Xp  —  0  (p  — 
mm{M,N}).  The  matrices  U  (M  x  M)  and  V(iV  x  N) 
are  formed  from  the  columns  and  {v,}^.j,  which 

are  the  orthogonal  eigenvectors  of  GGT  and  GTG,  respec¬ 
tively  [13].  Using  (7),  the  MNLS  solution  to  (5)  is  then 
given  by 


.  V^f 

HMNLS  =  >  - Vi 

i= 1  * 


(8) 


Note  that,  if  G  is  full  rank,  the  MNLS  is  equivalent  to  the 
LS  solution  in  the  overdetermined  case  and  to  the  MN  solu¬ 
tion  in  the  underdetermined  case.  For  example,  for  M  <  N 
(underdetermined),  1imnls  =  hAW  =  GT(GGr)_1f. 


where  ||  •  ||2  is  the  h  norm. 

In  the  context  of  the  linear  and  quadratic  prediction  ap¬ 
proach  taken  in  this  paper,  the  MSE  decreases  monotoni- 
cally  with  k.  A  criterion  for  choosing  an  appropriate  value 
of  k  is  needed.  In  the  context  of  contrast-agent  imaging, 
an  obvious  criterion  is  the  contrast.-to-tissue  ratio  (CTR) 

CTR  =  lOloglo0Q,  (11) 

where  Pc  is  the  average  power  of  signals  in  a  UCA  region, 
and  Pt  is  the  average  power  of  signals  in  a  tissue  region. 
The  average  power  of  signals  in  a  given  region,  P,  can  be 
expressed  as 

P  =  TJ  '  (12) 

where  Xij  is  the  signal  in  that  region.  Note  that  the  CTR 
used  in  this  algorithm  is  determined  from  quadratic  com¬ 
ponents  of  the  SVF  model.  The  CTR  provides  a  meaning¬ 
ful  stopping  criterion  for  TSVD  since  E(k)  is  monotoni- 
cally  decreasing  or  nonincreasing  with  k. 

One  may  use  a  variety  of  methods  to  obtain  a  regular¬ 
ized  solution  to  the  set  of  equations  involving  the  Volterra- 
kernels.  Examples  are  linear  or  nonlinear  programming  in 
cases  when  certain  constraints  on  the  filter  coefficients  are 
known  to  apply  (e.g.,  positivity)  or  penalized  maximum 
likelihood  when  known  statistical  properties  of  the  data 
can  be  incorporated.  A  powerful  approach  seeks  a  solution 
to  a  constrained  optimization  problem  of  the  form 

min  B?  subject  to  Gh  =  f,  (13) 

hq 
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where  R 2  is  an  appropriately  chosen  quadratic  ratio  to  be 
minimized  (e.g.,  inverse  of  CTR).  For  example,  R 2  can  be 
chosen  to  reflect  the  inverse  of  the  CTR  before  log  com¬ 
pression.  In  general,  it  can  be  any  quantity  that  depends 
on  the  solution  hq  that  is  to  be  minimized.  This  leads  to 
a  modified  form  of  the  TSVD  given  here 


+  1-R? 


uffVi 


(14) 


where  7  is  an  appropriately  chosen  threshold,  and  R'f  is 
quadratic,  ratio  resulting  from  the  quadratic,  kernel  ob¬ 
tained  from  the  ith  singular  mode.  Regardless  of  the  reg¬ 
ularization  procedure,  the  fundamental  result  here  is  the 
use  of  the  linear/quadratic  prediction  to  obtain  a  set  of  in¬ 
dependent  equations  that  can  be  solved  robustly  to  obtain 
the  Volterra  filter  kernels. 


C.  Quadratic  Images 


Quadratic  images  are  obtained  from  quadratic  compo¬ 
nents  of  the  second-order  Volterra  model.  The  coefficients 
of  the  SVF  are  derived  from  the  beamformed  RF  data 
taken  from  a  representative  region  on  a  standard  B-mode 
image.  Details  of  the  algorithm  to  produce  the  quadratic 
image  are  as  follows. 

From  the  standard  B-mode  image,  a  UCA  region  and 
a  tissue  region  are  defined  for  the  CTR  computation.  The 
definition  of  the  CTR  reference  regions  depends  on  the 
imaging  target  are  described  in  Section  III-C.  In  general, 
we  try  to  find  regions  at  the  same  depth  and  with  the 
same  beam  angle  with  respect  to  the  axis  of  the  imaging 
array.  Furthermore,  whenever  possible,  we  chose  multiple 
overlapping  subregions  to  obtain  multiple  CTR  values  at 
different  depths. 

Once  the  CTR  reference  regions  are  defined,  a  segment 
of  RF  data  from  an  axial  line  is  selected  to  form  a  sys¬ 
tem  of  linear  equations  according  to  (4).  This  segment 
can  be  selected  from  the  tissue  or  the  UCA  region  as  long 
as  the  appropriate  regularization  of  the  MNLS  solution  is 
sought.  For  example,  when  TSVD  is  used  for  regulariza¬ 
tion,  CTRs  of  quadratic  signals  calculated  from  various  or¬ 
ders  of  TSVD  solutions  are  collected.  With  a  defined  range 
of  system  orders,  a  CTR  plane  as  a  function  of  truncation 
parameters  and  system  orders  is  determined.  Filter  coef¬ 
ficients  for  the  quadratic  imaging  generation  are  obtained 
from  a  truncation  parameter  and  a  system  order  that  give 
the  highest  CTR  value  in  the  CTR  plane.  Of  course,  the 
TSVD  approach  can  be  used  to  obtain  the  coefficients  of 
the  quadratic  kernel,  h,Q{i,j),  for  a  predetermined  filter 
order  m.  This  may  be  necessary  from  the  implementation 
point  of  view,  when  the  size  of  the  kernel  is  to  be  kept  at 
manageable  level.  The  results  presented  in  this  paper  are 
obtained  with  low-order  filter  to  emphasize  the  practicality 
of  the  Volterra  filter  approach. 

The  quadratic  image  is  produced  by  applying  the 
quadratic  filter  coefficients  to  the  beamformed  RF  data 


Fig.  2.  A  flowchart,  of  the  algorithm  for  quadratic  image  generation. 


throughout  the  standard  B-mode  image  to  estimate  the 
quadratic  component 


777  — 1  in—  1 

UQ(n  +  1)  =  S  u(n  ~  -  fc)- 

i— 0  k=j  (15) 

where  tiQ(j ,  k )  is  the  estimated  quadratic  kernel  (extracted 
from  Iimnls)-  A  flowchart  of  this  algorithm  is  shown  in 
Fig.  2. 

D.  Harmonic  Imaging 

Second  harmonic  imaging  is  based  on  filtering  the  RF 
data  with  a  zero-phase  linear  BPF  centered  at  twice  the 
fundamental  with  restricted  bandwidth  to  minimize  the 
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overlap  between  the  fundamental  and  the  second  har¬ 
monic.  The  approach  is  appropriate  for  native  harmonic 
imaging.  However,  when  imaging  UCA,  the  echo  data 
from  the  UCA  regions  tends  to  have  wider  bandwidth 
when  compared  with  echoes  from  tissue  regions.  There¬ 
fore,  to  obtain  the  highest  possible  contrast  with  a  lin¬ 
ear  harmonic  filter,  we  varied  both  the  center  frequency 
and  the  bandwidth  of  the  filter  and  computed  a  CTR 
plane  as  a  function  of  these  parameters.  This  approach 
produced  superior  results  compared  to  standard  SH  imag¬ 
ing.  Harmonic  imaging  results  shown  in  Section  IV  are 
obtained  with  the  optimal  linear  BPF  designed  using  the 
Parks-McClellan  algorithm  [18].  All  linear  filters  used  in 
this  paper  were  implemented  in  a  zero-phase  realization 
(y(n)  =  h(n)  *  h(—n)  *  x(n)  or  Y{f)  =  \H(f)\2X(f)).  The 
order  of  the  filter  was  chosen  to  achieve  50  dB  stopband 
attenuation  and  0.5  passband  ripple. 


III.  Materials  and  Methods 

Imaging  results  shown  below  demonstrate  the  contrast 
enhancement  achieved  by  the  postbeamforming  SVF  in  the 
context  of  contrast-agent  imaging  in  tissue-mimicking  me¬ 
dia.  Two  imaging  targets  with  two  imaging  probes  are  used 
to  demonstrate  generality  and  robustness  of  the  approach. 

A.  Contrast  Agent 

The  contrast  agent,  BR14  (Bracco  Research  S.A., 
Geneva,  Switzerland),  was  used.  The  BR14  is  a  new  ex¬ 
perimental  agent,  that  consists  of  high  molecular  weight 
gas  bubbles  encapsulated  by  a  flexible  phospholipid  shell. 
A  0.125  mL  sample  of  BR14,  prepared  with  5  mL  of  0.9% 
saline,  was  diluted  in  500  mL  of  0.9%  saline  leading  to 
a  1:4000  dilution.  In  addition  to  BR14,  cellulose  particles 
(Sigma  Cell  Type  20,  Sigma  Chemical  Co.,  St.  Louis,  MO) 
were  used  as  linear  scatterers  in  the  flow  channel  target  de¬ 
scribed  below. 

B.  Imaging  Targets 

Images  of  two  targets  from  two  experimental  setups 
were  used  to  evaluate  the  performance  of  the  quadratic 
imaging  technique.  The  first  experimental  setup  is  shown 
in  Fig.  3.  The  target  is  the  L-shaped  tissue-mimicking 
phantom  in  a  beaker  containing  a  dilution  of  BR14. 
The  BR14  was  constantly  stirred  by  a  magnetic  stir¬ 
rer.  A  MEGAS  scanner  (ESAOTE  S.p.A.,  Genova,  Italy) 
equipped  with  a  2-MHz  phased  array  probe  was  used  for 
image  acquisition.  The  RF  data  acquisition  was  performed 
with  15-bit  resolution  at  40-MHz  sampling  frequency  and 
without  time  gain  control  (TGC).  This  experiment  allowed 
us  to  compare  echoes  from  three  different  locations:  the 
L-shaped  tissue-mimicking  phantom,  the  BR14,  and  the 
echo-free  region  visible  in  the  scan.  The  significance  of  the 
latter  is  that  any  signal  components  observed  in  this  region 


Tissue  M'imickingl 
[L-Shaped  Phantoni 


Beaker  [■ 


1 


Phased  Array 


BR14 
Contrast  Agent! 


Echo-Free  Region 
Visible  in  the  Scan 


Fig.  3.  Schematic  for  the  imaging  setup  for  the  L-shaped  phantom 
surrounded  by  the  UCA. 


Fig.  4.  The  imaging  setup  for  the  flow  phantom. 


are  largely  artifacts  from  beamforming  and/or  reverbera¬ 
tion.  It  should  be  noted,  however,  that  there  may  be  a 
small  backscatter  component  present  due  to  the  beam  re¬ 
flections  at  the  beaker,  but  this  is  hard  to  quantify  in  the 
presence  of  the  contrast  agent.  Despite  this  limitation,  it 
is  interesting  to  compare  the  nature  of  the  image  pixels  in 
this  region  from  the  three  imaging  methods. 

The  second  experimental  setup,  shown  in  Fig.  4,  was 
used  in  obtaining  images  of  a  flow  phantom  (Model  524; 
ATS  Laboratories,  Inc.,  Bridgeport,  CT)  containing  four 
flow  channels  with  diameters  2,  4,  6,  and  8  mm  embedded 
in  rubber-based  tissue  mimicking  material.  The  flow  phan¬ 
tom  was  connected  to  a  flow  system  with  a.  roller  pump 
(Model  77200-60;  Cole-Parmer  Instrument  Co.,  Vernon 
Hills,  IL).  Subsequently,  the  diluted  BR14  and  cellulose 
particles  were  circulated.  In  addition,  the  diluted  BR14 
was  constantly  stirred  in  a  beaker  using  a  magnetic  hot 
plate  stirrer  (EW-84303-20;  Corning  Inc.,  Corning,  NY). 

This  experiment  was  designed  to  compare  linear 
backscattered  signals  from  cellulose  in  the  8-mm  diame¬ 
ter  flow  channel  with  nonlinear  backscattered  signals  from 
BR14  in  the  6-mrn  diameter  flow  channel.  The  RF  data 
were  recorded  and  saved  for  later  processing  by  the  Tech- 
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nos  MP  ultrasound  system  (ESAOTE  S.p.A.,  Genova, 
Italy)  with  a  convex  array  probe  (CA421;  ESAOTE  S.p.A., 
Genova,  Italy)  located  perpendicularly  to  both  the  UCA 
and  the  cellulose  flow  channels.  The  RF  data  were  acquired 
with  16-bit  resolution  at  20-MHz  sampling  frequency  and 
without  TGC.  A  3-cycle  pulse  at  3.13  MHz  was  transmit¬ 
ted  to  produce  standard  B-mode  images. 


C.  Quantitative  Analysis 

After  the  RF  data  corresponding  to  standard  B-mode 
images  were  collected  from  experimental  setups  described 
above,  quadratic  images  were  obtained  using  the  algorithm 
described  in  Section  II-C.  For  comparison,  harmonic  im¬ 
ages  were  produced  by  filtering  the  standard  B-mode  im¬ 
age  with  bandpass  filters  with  center  frequencies  and  frac¬ 
tional  bandwidths  chosen  to  maximize  the  CTR  (defined 
by  (11))  as  described  in  Section  II-D. 


B-mode  and  the  output  of  the  harmonic  or  the  quadratic 
filter  for  the  harmonic  and  quadratic  images,  respectively. 

Amin  and  Amax  are  minimum  and  maximum  values,  re¬ 
spectively.  Both  Ami„  and  Amax  are  carefully  computed 
so  that  histograms  of  images  cover  256  gray  levels  without 
saturation.  This  was  achieved  by  finding  the  average  of  the 
minima  (maxima)  of  all  image  lines  after  median  filtering. 
This  was  done  to  ensure  that  the  display  range  is  not  set 
by  extreme  values  of  Amjn  or  Amax. 

Rom  gray-level  images,  the  contrast  ratio  (CR)  [15] 
used  as  the  contrast  measurement  between  any  two  regions 
is  given  by 


CR  = 


\h  ~  h\ 
V*  l+e'f 


(18) 


where  and  I2  are  the  average  of  gray  levels,  and  o\  and 
o-i  are  the  corresponding  standard  deviations  in  the  first 
region  and  the  second  region,  respectively. 


1.  CTR  Computations:  The  CTR  was  used  for  com¬ 
parison  between  the  image  data  in  the  RF  domain  before 
scan  conversion  as  described  in  (11).  Size  and  location  of 
reference  regions  for  each  target  are  as  follows.  For  the 
L-shape  phantom,  45-mrn  axial  segments  of  20  adjacent 
A-lines  on  the  tissue  side  and  the  same  number  on  the 
contrast  side  were  used.  This  phantom  provided  an  excel¬ 
lent  opportunity  to  compute  CTR  at  multiple  regions  in 
the  contrast  and  tissue  regions  with  the  same  depth  and 
the  same  beam  angle  from  the  axis  of  the  imaging  array. 
The  calculation  region  extends  by  45  mm  in  the  axial  di¬ 
rection  (from  50  mm  to  95  mm)  and  comprises  20  adjacent 
A-lines  of  RF  data  before  scan  conversion.  The  CTR.  val¬ 
ues  are  calculated  from  cells  with  3.8-mm  axial  extent  with 
10  connected  A-lines  and  50%  overlap  in  both  directions. 
For  the  flow  phantom,  one  contrast  region  is  in  UCA  flow 
channel  and  the  tissue  region  located  between  two  flow 
channels  with  the  same  axial  extent  (3  mm  in  this  case). 
Twenty  A-line  segments  from  the  each  region  were  used. 

2.  Contrast  Ratio:  As  can  be  seen  from  results  given 
below,  the  dynamic  range  (in  decibels)  of  images  from 
quadratic  components  is  approximately  twice  the  dynamic 
range  of  standard  B-mode  and  SH  images.  In  order  to 
account  for  image  perception  on  standard  8-bit  display, 
all  images  are  represented  with  their  full  dynamic  range 
mapped  to  256  gray  levels.  That  is, 


3.  Histogram,  and  Receiver  Operating  Characteristic 
Analysis:  In  the  imaging  results  shown  below,  we  evaluate 
the  CTR  and  CR  based  on  regions  in  the  image  represen¬ 
tative  of  the  UCA  and  tissue  regions  with  the  same  num¬ 
ber  of  pixels  and  at  the  same  depth.  When  the  number  of 
pixels  in  these  regions  is  sufficiently  high  to  produce  mean¬ 
ingful  statistics,  histograms  of  the  pixel  data  are  produced 
to  demonstrate  these  statistics.  In  addition,  receiver  oper¬ 
ating  characteristics  (ROC)  analysis  [16]  is  performed  as  a 
simple  method  for  classification  of  the  different  regions.  In 
addition,  when  additional  regions  can  be  identified  (e.g., 
low-scattering  regions)  in  any  image,  the  CR  between  tis¬ 
sue  and  such  regions  is  also  evaluated.  Please  note  that  the 
ROC  curves  are  included  in  this  paper  to  further  quantify 
the  degree  of  overlap  between  different  histograms  from 
tissue,  contrast,  and  low-scattering  regions. 


D.  Resolution  Measurements 

We  used  correlation  lengths  calculated  from  the  2D  au¬ 
tocorrelation  of  echoes  in  tissue  regions  as  described  in  [17] 
to  measure  spatial  resolution.  Intensity  images  are  scan 
converted  and  uniform  speckle  regions  were  identified  to 
compute  the  average  speckle  correlation  cell  size 


1  = 


A  -  A, 


iin  \ 
hn in  / 


x  255, 


(16) 


= 


Ci(x,y ) 
<5/(0, 0) 


dx  dy , 


(19) 


where  A  is  the  log  magnitude  pixel  values,  given  by 


A  =  20  log 


/  \n(x)}  ) 

10  I  [  ^ref  J  J  ’ 


(17) 


where  Tt(x)  denotes  the  Hilbert  transform  of  data  vector, 
2,  and  a;ref(>  0)  is  some  constant  (typically  maximum  am¬ 
plitude).  The  data  vector  is  the  digitized  RF  for  standard 


where  Ci(x,y)  is  the  2D  correlation  function  of  the  inten¬ 
sity  autocovariance  function,  and  X  and  Y  are  taken  to 
be  sufficiently  large  to  allow  the  magnitude  of  the  auto¬ 
covariance  to  drop  to  negligible  levels.  The  vertical  and 
horizontal  correlation  cell  sizes,  representing  the  axial  and 
lateral  resolution,  respectively,  were  used  to  compare  spa¬ 
tial  resolution  for  the  three  imaging  techniques. 
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Fig.  5.  Average  spectra  from  the  contrast  (thick)  and  tissue  (thin) 
regions  of  the  L-shaped  phantom. 

IV.  Results 

A.  L-Shaped  Phantom 

Fig.  5  shows  the  average  spectra  of  typical  echoes  from 
the  contrast  and  tissue  regions  of  the  L-shaped  phantom. 
The  average  spectra  were  calculated  by  averaging  win¬ 
dowed  periodogram  of  every  echo  line  in  regions  described 
in  Section  III-C.l.  These  spectra  suggest  that  echoes  from 
the  contrast  region  have  broader  bandwidth  compared  to 
the  echoes  from  tissue-like  region.  As  mentioned  in  Sec¬ 
tion  II-D,  this  suggests  that  the  best  linear  filter  for  con¬ 
trast  enhancement  may  not  be  the  standard  SH  filter. 
Rather,  a  general  BPF  with  center  frequency  and  relative 
bandwidth  designed  to  maximize  the  CTR  may  be  sought. 
Based  on  this  approach,  we  have  designed  a  zero-phase 
BPF  with  center  frequency  of  3.2  MHz  and  relative  band¬ 
width  of  30%  to  demonstrate  the  performance  of  linear 
postbeamforming  harmonic  imaging.  The  quadratic  filter 
was  derived  from  echoes  from  the  tissue  region  according  to 
the  algorithm  described  in  Section  II-C  (with  P  =  Q  =  10 
in  Fig.  2).  Panels  (a),  (b),  and  (c)  of  Fig.  6  show  images  ob¬ 
tained  using  standard  B-mode,  linear  BPF,  and  quadratic 
imaging  methods,  respectively.  Due  to  differences  in  dy¬ 
namic  ranges  for  the  three  methods,  each  image  is  dis¬ 
played  with  its  full  dynamic  range  as  can  be  seen  from 
the  decibel-level  scale  bars  shown  in  accordance  with  (16). 
Fig.  6(d)  shows  four  rectangular  regions  used  for  charac¬ 
terization  of  the  imaging  results.  Regions  A\  and  are 
representative  of  tissue  and  contrast  echoes,  respectively. 
Regions  B\  and  B2  are  representative  of  tissue  and  echo- 
free  region,  respectively.  These  regions  are  used  for  CTR 
calculations  as  well  as  evaluating  histograms  of  the  recon¬ 
structed  images.  Due  to  the  structure  of  the  target,  regions 
B\  and  B2  could  not  be  placed  at  the  same  depth  in  the  im¬ 
age.  Nevertheless,  region  B\  is  representative  of  tissue  re¬ 


sponse  and  provides  a  valid  baseline  for  comparison  of  con¬ 
trast  with  the  echo-free  region,  B2 •  The  standard  B-mode 
image  (Fig.  6(a))  was  obtained  from  the  digitized  RF  echo 
data  and  displayed  with  a  dynamic,  range  of  60  dB.  The 
structure  of  the  L-shaped  phantom  can  be  recognized  with 
low  contrast  between  the  UCA  and  the  tissue-mimicking 
regions.  The  strong  specular  reflection  at  the  top  of  the 
UCA  region  is  due  to  a  boundary  layer  formed  by  the 
agent  at  the  interface  with  the  tissue-mimicking  medium. 
A  careful  examination  of  the  image  reveals  a  higher  level 
of  contrast  between  the  two  regions  at  close  range.  A  CTR 
value  of  6.65  dB  was  computed  based  on  echo  signals  from 
regions  A2  and  A\  (Fig.  6(d)). 

The  linear  harmonic  BPF  image  shown  in  Fig.  6(b)  was 
produced  by  filtering  the  RF  data  with  bandpass  filter  cen¬ 
tered  around  the  second  harmonic  (3.2  MHz)  with  a  frac¬ 
tional  bandwidth  of  30%.  The  center  frequency  and  the 
fractional  bandwidth  of  the  harmonic  filter  were  chosen  to 
optimize  the  CTR  value  for  regions  A\  and  A 2  in  Fig.  6(d). 
The  CTR  value  for  the  harmonic  image  (18.2  dB)  is  con¬ 
sistent  with  the  perceived  enhancement  clearly  visible  in 
the  image.  One  can  also  observe  the  loss  of  the  specular 
reflections  at  the  top  right  and  bottom  left  of  the  harmonic, 
image.  This  is  typical  since  these  echoes  have  significant 
low  frequency  components.  As  expected,  the  speckle  in  the 
tissue  region  appears  finer  than  that  of  the  standard  B- 
mode. 

The  quadratic  image  obtained  using  the  algorithm  de¬ 
scribed  in  Section  II-C  is  shown  in  Fig.  6(c).  The  CTR 
for  this  image  is  21.3  dB  indicating  contrast  enhancement 
over  both  standard  B-mode  and  harmonic  images.  One 
important  feature  of  the  quadratic  image  is  the  increased 
dynamic  range  compared  with  B-mode  and  harmonic  im¬ 
ages.  This  increased  dynamic  range  results  in  contrast  en¬ 
hancement  without  loss  in  image  features,  such  as  specular 
reflections,  which  may  be  of  diagnostic  value,  in  some  cases. 

Another  interesting  comparison  among  the  three  im¬ 
age  results  shown  in  Fig.  6(a.)-(c.)  is  the  echo-free  region 
visible  at  the  right  edge  of  the  scan.  One  can  see  that, 
contrast  between  tissue  and  this  region  in  the  quadratic 
image  is  superior  to  that  from  the  B-mode  image,  whereas 
this  contrast  in  the  harmonic  image  is  the  lowest.  This 
result  is  significant  because  echo  signals  from  this  region 
are  largely  artifacts  (due  to  beamforming  and/or  reverber¬ 
ations).  Quantitative,  measures  of  this  contrast  enhance¬ 
ment  between  the  tissue  and  echo-free  regions  are  given 
below. 

To  further  illustrate  the  imaging  results  given  in  Fig.  6, 
vertical  and  horizontal  lines  through  images  a, re  plotted 
in  Fig.  7.  Axial  lines  from  the  three  imaging  techniques 
through  the  UCA  region  are  shown  in  Fig.  7(a).  One  can 
see  that  the  specular  reflector  is  all  but  eliminated  in  the 
harmonic  image,  while  it  remains  visible  for  the  other  two 
methods.  Further,  the  axial  line  from  the  quadratic  image 
shows  no  apparent  loss  in  the  axial  resolution  at  the  spec¬ 
ular  reflection  location.  Lateral  lines  at  the  90- mm  depth 
from  the  three  imaging  techniques  are  shown  in  Fig.  7(b). 
These  lateral  lines  pass  through  the  UCA,  tissue,  and  the 
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Fig.  6.  Images  of  the  beaker  containing  the  tissue-mimicking  L-shaped  phantom  surrounded  by  the  UCA  from  three  imaging  techniques, 
(a)  Standard  B-mode.  (b)  Linear  harmonic  BPF.  (c)  Quadratic,  (d)  Boxes  indicate  regions  for  CTR  and  histogram  calculations. 


echo-free  regions.  One  can  see  the  contrast  enhancement 
from  the  lateral  line  of  the  harmonic  image  compared  with 
the  lateral  line  of  the  standard  B-mode  image.  Most  signif¬ 
icantly,  however,  one  can  see  that  the  quadratic  image  pro¬ 
vides  the  maximum  separation  between  the  three  regions 
without  apparent  loss  of  lateral  resolution.  This  was  con¬ 
firmed  by  computing  the  correlation  cell  size  in  the  scan- 
converted  intensity  images  according  to  (19)  [17]  from  the 
three  imaging  methods  in  the  tissue  region  (59  to  73  mm 
axial  and  -20  to  —14  mm  lateral  in  Fig.  6).  These  values 
are  reported  in  Table  I. 

B.  Flow  Channel  Phantom , 

Fig.  8  shows  the  average  spectra  of  typical  echoes  from 
the  contrast  and  tissue  regions  of  the  flow  phantom  de¬ 
scribed  in  Section  III-B.  The  average  spectra  are  calcu¬ 
lated  by  averaging  windowed  periodogram  of  every  echo 


line  in  regions  described  in  Section  III-C.l.  As  with  the 
L-shaped  phantom  result,  the  echoes  from  the  UCA  re¬ 
gion  exhibit  broader  bandwidth  than  those  from  tissue  re¬ 
gion.  The  standard  B-mode  image  of  the  flow  phantom 
consisting  of  the  UCA  and  cellulose  in  flow  channels  was 
acquired  using  the  experimental  setup  described  in  Sec¬ 
tion  III-B.  The  image  is  shown  in  Fig.  9(a).  One  can  see 
the  backscattered  enhancement  due  to  the  UCA  in  the  6- 
nnn  flow  channel  (55-  to  61-mm  axial  and  —6  to  0  mm 
lateral)  compared  with  the  surrounding  tissue  with  CTR. 
9.8  dB.  The  CTR  was  calculated  based  on  echo  signals 
from  the  6-mm  flow  channel  excluding  the  specular  reflec¬ 
tors  (Box  Ci)  and  echoes  from  tissue  at  the  same  depth 
(Box  C2T  On  the  other  hand,  the  echoes  from  the  cellulose 
in  the  8-mm  flow  channel  (53  to  61  mm  axial  and  17  to 
25  mm  lateral)  are  weaker  than  those  from  the  surrounding 
tissue.  A  linear  scatterer-to-tissue  ratio  (LTR.)  was  calcu¬ 
lated  based  on  echo  signals  from  the  8-mm  flow  channel 


PHUKPATTARANONT  AND  EBBINI:  POST-BEAMFORMING  SECOND-ORDER  VOLTERRA  FILTER 


995 


TABLE  I 

CR,  CTR,  AND  Correlation  Cell  Size  Values  for  the  L-Shaped  Imaging  Target.  The  CTR  Values  are  Given  as  Mean  in 

Decibels  with  Standard  Deviation  in  Parentheses. 


Imaging  method 

UCA/T-Mimic 

T-Mimie./E-F 

CTR 

S* 

Sx 

B-mode 

0.94 

1.51 

6.65  (1.3) 

1.2 

2.3 

SH 

2.56 

1.12 

18.2  (1.6) 

1.0 

1.6 

Quadratic 

2.16 

1.66 

21.3  (1.9) 

0.42 

0.8 

(a) 


(b) 

Fig.  7.  Lines  through  images  in  Fig.  6  from  three  imaging  techniques. 
Thin:  B-mode.  Dash:  Harmonic.  Thick:  Quadratic.  (Top)  Axial  lines 
through  the  UCA  region.  (Bottom)  Lateral  lines  (at  depth  90  mm) 
through  the  tissue-mimicking,  UCA,  and  echo-free  regions. 


Fig.  8.  Average  spectra  from  the  contrast  (thick)  and  tissue  (thin) 
regions  of  the  flow  phantom. 


(Box  C3)  and  echoes  from  tissue  at  the  same  depth  (Box 
C2).  The  LTR  was  -8.8  dB,  calculated  based  on  (11)  for 
8-mrn  channel. 

The  harmonic,  image  shown  in  Fig.  9(b)  was  produced 
by  filtering  the  RF  data,  with  bandpass  filter  centered 
around  4.8  MHz  with  a  fractional  bandwidth  of  30%.  The 
center  frequency  and  the  fractional  bandwidth  of  the  har¬ 
monic  filter  were  chosen  to  optimize  the  CTR  value  using 
the  method  described  in  Section  II-D  for  the  UCA  and 
tissue  regions  shown  in  Fig.  9(d).  The  CTR.  value  for  the 
harmonic  image  is  15  dB,  which  is  consistent  with  the  per¬ 
ceived  enhancement  clearly  visible  in  the  displayed  image. 
It  is  interesting  to  note  the  loss  of  the  specular  reflections 
from  the  contrast  flow  channel  but  not  from  the  linear  scat- 
terer  channel.  It  is  also  interesting  to  note  that  the  speckle 
in  the  tissue  region  is  finer  than  that  of  the  standard  B- 
mode.  This  indicates  that  the  optimization  of  the  harmonic 
filter  resulted  in  a  larger  fractional  bandwidth  than  that 
of  the  fundamental  component.  The  LTR  for  the  harmonic 
component  is  now  —7.6  dB.  That  is,  the  harmonic  image 
does  not  preserve  the  low-scattering  region  as  it  improves 
the  contrast  between  the  UCA  and  tissue  regions. 

Fig.  9(c)  shows  the  quadratic  image  obtained  using  the 
algorithm  described  in  Section  II-C  (with  P  =  Q  —  12 
in  Fig.  2).  Compared  with  the  standard  B-mode,  the  en¬ 
hanced  visualization  of  the  UCA  in  the  small  flow  channel 
surrounded  by  the  tissue  can  be  seen.  This  enhanced  visu¬ 
alization  is  also  consistent  with  the  higher  CTR  (22.2  dB). 
In  addition,  one  can  see  the  reduction  in  the  echogenic¬ 
ity  of  the  cellulose  in  the  8-mm  flow  channel.  This  is  also 
reflected  in  the  value  of  LTR  (—15.4  dB)  computed  for 
this  channel.  This  result  demonstrates  the  ability  of  the 
quadratic  filter  to  separate  the  tissue  from  both  the  UCA 
and  the  low-scattering  cellulose  region  simultaneously.  We 
can  also  see  the  preservation  of  specular  reflections  from 
both  the  UCA  and  cellulose  flow  channels.  This  can  be  seen 
more  quantitatively  in  Fig.  10,  which  shows  the  decibel 
values  for  axial  and  lateral  lines  through  images  in  Fig.  9. 
Axial  lines  through  the  center  of  the  cellulose  and  the  UCA 
flow  channels  are  shown  in  Fig.  10(a.)  and  Fig.  10(b),  re¬ 
spectively.  Lateral  lines  through  the  center  of  both  the  cel¬ 
lulose  and  the  UCA  flow  channels  are  shown  in  Fig.  10(c). 
One  can  see  the  contrast  enhancement  in  different  me¬ 
dia  (i.e. ,  the  UCA,  the  cellulose,  and  the  tissue)  from  the 
quadratic  image  without  loss  in  spatial  resolution.  In  par¬ 
ticular,  for  the  B-mode  image,  the  signal  from  the  contrast 
channel  is  13  dB  above  the  tissue  level  and  30  dB  above 
the  cellulose  channel.  On  the  other  hand,  for  the  quadratic 
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Fig.  9.  Images  of  the  flow  phantom  using  (a)  Standard  B-mode.  (b)  Harmonic,  (c)  Quadratic,  (d)  Boxes  Ci.Ca,  and  C3  indicate  regions  for 
CTR  and  CR  calculations.  Note  the  hyperechoic  region  from  the  UCA  in  the  6-mm  flow  channel  (55  mm  to  61  mm  axial  and  —6  mm  to  0 
lateral)  and  the  hypoechoic  region  from  the  cellulose  (53  mm  to  61  mm  axial  and  17  mm  to  25  mm  lateral). 


image,  the  signal  from  the  contrast  channel  is  30  dB  higher 
than  tissue  and  67  dB  higher  than  the  cellulose  channel.  Fi¬ 
nally,  examination  of  the  strong  specular  reflections  along 
the  axial  lines  shown  in  Fig.  10  suggests  that  the  quadratic 
filter  preserves  the  resolution  of  the  system.  This  was  con¬ 
firmed  by  computing  the  correlation  cell  size  in  the  scan- 
converted  intensity  images  according  to  (19)  [17]  from  the 
three  imaging  methods  in  the  tissue  region  (52  to  65  mm 
axial  and  4  to  17  mm  lateral  in  Fig.  9).  These  values  are 
reported  in  Table  II. 

C.  Classification  Results 

To  give  the  reader  a  quantitative  idea  of  the  ability  of 
the  three  different  imaging  methods  to  separate  the  dif¬ 
ferent  regions  in  the  imaging  targets,  we  use  histograms 
from  representative  areas  of  these  regions.  For  example, 


for  the  tissue-mimicking  L-shaped  phantom  surrounded 
by  the  UCA  in  the  beaker,  regions  from  different  media 
are  defined  in  Fig.  6(d)  as  follows:  the  region  A1  and  A2 
(the  tissue  and  UCA  regions,  respectively)  and  the  re¬ 
gion  B1  and  B2  (the  tissue  and  echo-free  regions,  respec¬ 
tively).  After  images  from  three  different  imaging  tech¬ 
niques  are  scan  converted  and  represented  with  256  levels 
of  gray  covering  their  full  dynamic  range,  the  histogram 
of  each  region  is  determined  and  shown  in  Fig.  11.  His¬ 
tograms  from  regions  A1  and  A2  are  shown  in  Fig.  11 
(left-hand  side).  One  can  see  the.  degree  of  overlap  be¬ 
tween  the  histograms  is  highest  for  the  standard  B-mode 
image.  On  the  other  hand,  the  corresponding  histograms 
are  well  separated  for  the  harmonic  and  quadratic  images. 
This  is  consistent  with  the  increased  level  of  contrast  per¬ 
ceived  from  direct  visualization  of  the  images  in  Fig.  6. 
Histograms  produced  from  regions  B1  and  B2  are  shown  in 
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Fig.  10.  Lines  through  the  standard  B-mode  (thin),  harmonic  (dashed),  and  the  quadratic  (thick)  from  images  shown  in  Fig.  9.  (Top)  Axial 
linos  through  the  center  of  the  cellulose  flow  channel.  (Middle)  Axial  lines  through  the  center  of  the  UCA  flow  channel.  (Bottom)  Lateral 
lines  through  the  center  of  both  two  flow  channels. 


TABLE  II 

CR,  CTR,  and  Correlation  Cell  Size  Values  for  the  Flow  Channel  Target.  The  CTR  Values  are  Given  in  Decibels. 


Imaging  method 

UCA/Tissue 

Tissue/Cellulose 

CTR 

Sz 

Sx 

B-mode 

1.57 

1.35 

9.8 

1.0 

1.2 

SH 

2.12 

1.21 

15.0 

0.6 

0.9 

Quadratic 

2.46 

1.65 

22.2 

0.9 

1.1 

Fig.  11  (right-hand  side).  Gray-level  histograms  of  regions 
B1  and  B2  produced  from  harmonic  image  has  higher  de¬ 
gree  of  overlap  than  those  from  the  standard  B-mode  and 
the  quadratic  images.  This  is  further  quantified  by  ROC 
curves  [16]  obtained  from  histograms  of  region  A1  (Tis¬ 
sue)  and  A2  (UCA)  and  region  B1  (tissue)  and  B2  (echo- 
free)  are  shown  in  Fig.  12(a)  and  Fig.  12(b),  respectively. 
As  can  be  seen  in  Fig.  12(a),  the  Az  values  between  re¬ 
gion  A1  and  A2  from  the  standard  B-mode,  harmonic  and 
quadratic  images  are  0.8269,  0.9936,  and  0.9876,  respec¬ 
tively.  These  values  demonstrate  improved  classification 
performance  between  region  A1  and  region  A2  from  har¬ 
monic  and  quadratic  images  over  the  standard  B-mode  im¬ 
age.  However,  as  shown  in  Fig.  12(b),  the  Az  value  of  the 
quadratic  image  (0.9535)  presents  the  best  classification 
performance  of  tissue  and  echo-free  regions  among  three 
imaging  techniques,  while  the  classification  performance 
from  the  harmonic  image  ( Az  =  0.8724)  is  inferior  to  that 


from  the  standard  B-mode  image  {Az  =  0.9302).  These 
measurements  show  that  the  quadratic  image  provides  im¬ 
proved  separation  of  the  UCA  from  the  tissue  mimic  com¬ 
parable  with  SH  performance.  At  the  same  time,  signifi¬ 
cant  improvement  in  separation  between  tissue  mimic  and 
echo-free  regions  is  achievable  in  the  quadratic  image  over 
both  the  standard  B-mode  and  SH  images. 

Contrast  ratios  determined  from  corresponding  gray 
scale  images  in  Fig.  6  are  shown  in  Table  I.  For  contrast 
comparison  between  the  UCA  and  tissue  regions,  CRs  ob¬ 
tained  from  regions  A1  and  A2  demonstrate  the  contrast 
enhancement  of  the  harmonic  (2.56)  and  the  quadratic 
(2.16)  images  over  the  standard  B-mode  image  (0.94).  On 
the  other  hand,  one  can  see  that,  the  CR  between  tissue  and 
echo-free  regions  of  the  harmonic  image  (1.12)  is  inferior 
to  that  of  the  standard  B-mode  image  (1.51).  However,  the 
quadratic  image  gives  the  maximal  CR  (1.66).  These  CR 
values  agree  with  both  the  visualization  of  images  shown 
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Fig.  11.  Gray-level  histograms  produced  from  images  shown  in  Fig.  6.  Top:  B-mode.  Middle:  Harmonic.  Bottom:  Quadratic.  Histograms 
on  the  left-hand  side  are  produced  from  region  A1  and  A2.  Thick:  Region  A1  (Tissue).  Thin:  Region  A2  (Contrast).  Histograms  on  the 
right-hand  side  are  produced  from  region  B1  and  B2.  Thick:  Region  B1  (Tissue).  Thin:  Region  B2  (Air). 


Fig.  12.  The  ROC  curves  produced  from  gray-level  histograms  in  Fig.  11  of  three  imaging  techniques.  Thin:  B-mode.  Dash:  Harmonic.  Thick: 
Quadratic.  The  ROC  curve  obtained  from  the  UCA  and  tissue  regions  is  shown  on  the  left-hand  side.  The  ROC  curve  obtained  from  the 
echo-free  and  tissue  regions  is  shown  on  the  right-hand  side. 
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in  Fig.  6  and  the  degree  of  overlap  between  the  histograms 
shown  in  Fig.  11.  These  results  have  demonstrated  that 
the  quadratic  image  provides  the  perceived  contrast  en¬ 
hancement  not  only  between  the  UCA  and  tissue  regions 
but  also  between  the  tissue  and  the  echo-free  regions.  We 
note  here  that  similar  analysis  was  performed  for  the  flow 
channel  data  and  that  the  results  are  in  full  agreement 
with  the  results  shown  in  this  section.  The  CR,  CTR,  and 
correlation  cell  size  values  for  both  the  L-shaped  and  flow 
phantoms  are  summarized  in  Tables  I  and  II. 

Similarly,  images  of  the  UCA  and  cellulose  flow  channels 
in  the  flow  phantom  as  shown  in  Fig.  9  are  scaled  into  8-bit, 
gray-level  images.  Gray-level  histograms  between  the  UCA 
and  tissue  regions  and  between  the  cellulose  and  tissue 
regions  are  produced  and  the  corresponding  ROC  curves 
are  determined.  Although  the  number  of  data  from  each 
region  is  quite  small  for  meaningful  statistical  analysis, 
results  from  both  CRs  and  Az  show  tendency  that  the 
quadratic  image  enhances  the  contrast  and  provides  better 
classification  performance,  compared  with  the  standard  13- 
mode  image.  Table  II  gives  the  CR  values  resulting  from 
the  flow  channel  experiment. 

V.  Discussion 

We  have  introduced  a  nonlinear  postbeamforming  fil¬ 
tering  algorithm  for  ultrasonic  pulse-ec.ho  imaging  based 
on  the  Volterra  filter  model.  The  main  goal  of  this  paper 
was  to  introduce  the  mathematical  basis  for  deriving  the 
filter  coefficients  from  standard  beamformed  RF  data  ob¬ 
tained  by  commercial  scanners.  In  addition,  we  have  pre¬ 
sented  imaging  results  from  two  laboratory  contrast  tar¬ 
gets  to  illustrate  the  nature  of  the  gray  scale  images  ob¬ 
tained  with  the  quadratic  component  compared  with  stan¬ 
dard  echo  signals  and  harmonic  images.  Images  from  the 
quadratic  components  were  obtained  by  applying  a.  sin¬ 
gle  filter  derived  from  echoes  from  the  contrast  region.  We 
have  decided  to  do  this  to  demonstrate  that  the  method  is 
quite  robust  in  the  sense  that  the  derived  filter  is  applied 
throughout  the  image  to  produce  images  free  of  artificial 
inhomogeneity.  This  is  a  desirable  method  from  the  im¬ 
plementation  point  of  view,  especially  for  real-time  imple¬ 
mentation. 

For  a.  fair  comparison,  all  images  presented  in  this  paper 
were  normalized  to  their  full  dynamic  range  and  displayed 
using  256  levels  of  gray.  Classification  results  based  on 
histogram  characterization  of  contrast  and  tissue  regions 
show  that  the  quadratic  images  produce  nearly  twofold 
increase  in  CR  values  (compared  with  standard  echo  im¬ 
ages)  without  loss  of  image  features  (compared  to  har¬ 
monic  images).  Classifying  the  tissue  and  echo-free  regions 
demonstrated  the  improved  performance  with  respect  to 
harmonic  imaging.  This  is  probably  due  the  vulnerability 
of  harmonic  images  to  noise  and  beamforming  artifacts, 
especially  when  the  sidelobes  of  the  transmit  beam  are  in 
contrast  regions. 

Even  though  we  have  focused  on  the  comparison  be¬ 
tween  gray-scale  images,  there  are  some  interesting  prop¬ 


erties  of  the  quadratic  signal  components  that  may  reveal 
important  information  on  the  nature  of  the  objects  pro¬ 
ducing  the  echo  signals.  For  example,  the  frequency  com¬ 
ponent  at  f  in  the  quadratic  signal  component  is  a  result 
of  all  frequency  components  f\  and  f-2  from  the  echo  sig¬ 
nal  such  that  /i  +  f'2  =  /,  weighted  by  £Tq(/i,/2).  This 
frequency  coupling  results  in  quadratic  signal  components 
with  high  SNR  values  due  to  rejection  of  additive  noise. 
An  interesting  question  is  whether  the  quadratic  compo¬ 
nents  are  more  directly  related  to  the  composition  of  the 
echoes  in  terms  of  coherent  and  diffuse  scattering.  This 
may  be  quite  significant  in  improving  the  robustness  of 
motion  tracking  and  displacement  estimation  algorithms. 

In  this  paper,  the  filter  is  designed  based  on  nonlin¬ 
ear  echoes  from  the  tissue-like  medium.  However,  the  fil¬ 
ter  does  not  distinguish  between  nonlinearities  from  UCA 
and  those  from  tissue.  The  contrast  enhancement  is  mainly 
due  to  the  higher  level  of  quadratic  component  (relative  to 
the  RF  signal  level)  in  the  contrast  regions.  For  example, 
in  the  L-shaped  phantom,  on  average,  the  quadratic  sig¬ 
nal  is  30  dB  and  50  dB  below  the  RF  in  the  contrast  and 
tissue  regions,  respectively.  It  is  interesting  to  note  that, 
for  this  target,  the  quadratic  component  from  the  mag¬ 
netic  bead  at  the  bottom  of  the  image  is  only  3  dB  be¬ 
low  the  RF  signal.  This  suggests  that  the  quadratic  filter 
does  not  completely  reject  quadratic  echoes  from  specu¬ 
lar  reflectors,  even  though  it  may  still  discriminate  against 
them.  Therefore,  for  applications  in  which  the  quadratic, 
components  from  the  contrast  agents  are  weak  or  compa¬ 
rable  to  tissue  components,  a  method  for  separating  the 
two  types  of  nonlinearity  is  needed.  There  are  several  pos¬ 
sible  approaches  to  the  filter  design  problem  so  that  better 
separation  between  nonlinear  echoes  from  UCA  and  tissue: 

•  Adaptive  implementation  of  the  SVF  based  on  fast 
recursive  LS  approach  [12].  This  will  allow  for  the  op¬ 
timization  of  filter  coefficients  based  on  the  local  level 
of  nonlinearity. 

•  Higher  order  filters  (e.g.,  cubic)  that  may  be  more 
sensitive  to  UCA  nonlinearity  than  tissue  nonlinear¬ 
ity.  This  is  due  to  the  observation  that,  under  nor¬ 
mal  imaging  conditions,  tissue  nonlinearity  is  at  most 
quadratic. 

•  Synthesis  of  pulse  sequences  to  excite  contrast  mi¬ 
crobubbles  to  maximize  their  nonlinear  response.  This 
is  motivated  by  the  recent  trend  in  contrast-agent, 
imaging,  which  calls  for  the  use  of  super  low  values 
of  MI.  The  identified  quadratic  kernel  of  the  SVF  can 
be  used  for  the  design  of  these  waveforms. 

We  have  chosen  to  compare  quadratic  images  with  im¬ 
ages  obtained  using  linear  harmonic  filters  rather  than  sec¬ 
ond  harmonic  filters.  This  was  based  on  the  observation 
that  the  spectra  of  echoes  from  contrast  regions  are  typi¬ 
cally  broader  than  those  from  tissue  regions  as  can  be  seen 
from  Figs.  5  and  8.  The  application  of  a  strict  SH  filter 
based  on  the  bandwidth  from  tissue  component  typically 
produce  inferior  results  when  UCA  regions  are  present  in 
the  imaging  field.  To  illustrate  this  point,  Fig.  13  shows 
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Fig.  13.  LHS:  The  image  of  L-shaped  phantom  produced  using  BPF  centered  at  4  MHz  with  fractional  bandwidth  25%.  RHS:  The  image  of 
flow  channel  phantom  produced  using  BPF  centered  at  6  MHz  with  fractional  bandwidth  25%. 


SH  images  of  both  the  L-shaped  and  flow  phantoms  based 
on  the  spectra  from  tissue  regions.  These  images  are  vis¬ 
ibly  inferior  to  those  obtained  using  the  linear  harmonic 
approach  described  in  Section  II-D.  We  note  that  quanti¬ 
tative  analysis  of  these  images  confirms  the  inferiority  of 
the  SH  images  shown. 

The  filters  used  for  extracting  the  quadratic  compo¬ 
nents  for  the  two  data  sets  were  of  orders  10  and  12  for 
the  L-shaped  phantom  and  the  flow-channel  phantom,  re¬ 
spectively.  In  deriving  the  filters,  we  have  assumed  that 
hQ(j,k)  ~  liQ(k,j),  which  implies  that  the  filters  can  be 
implemented  with  65  and  90  independent  coefficients,  re¬ 
spectively.  On  the  other  hand,  the  linear  bandpass  filter 
had  77  coefficients.  This  number  is  effectively  doubled  by 
the  zero-phase  implementation  described  in  Section  II-D. 
While  it  is  not  our  objective  to  compare  the  computational 
efficiency  of  the  quadratic  filters  with  the  linear  filter,  we 
note  that  the  requirements  of  the  quadratic  filters  will  not 
present  a  severe  problem  for  real-time  implementation  on 
modern  ultrasound  scanners.  We  also  note  that  efficient 
software  and  hardware  implementations  of  the  Volterra  fil¬ 
ter  have  been  extensively  studied  in  recent  literature  [19]. 

Finally,  even  though  the  two  imaging  targets  in  this  pa¬ 
per  contain  UCA  regions,  quadratic  images  can  be  used  in 
“native,  quadratic”  form,  in  much  the  same  way  as  SH  im¬ 
ages.  The.  enhancement  of  the  contrast  between  the  tissue 
mimic  and  the  low  scattering  regions  in  both  targets  offer 
an  illustration  of  the  nature  of  quadratic  images  in  native 
mode.  We  are  currently  investigating  speckle  reduction  in 
quadratic  images  (some  speckle  reduction  can  be  observed 
in  the  tissue  mimic  in  quadratic  images  from  both  phan¬ 
toms).  This  is  the  subject  of  a  future  report. 

VI.  Conclusions 

A  nonlinear  postbeamforming  filter  based  on  the 
Volterra  model  was  introduced  in  this  paper.  The 


quadratic  component  from  the  SVF  signal  separation 
model  was  shown  to  produce  images  with  high  contrast 
and  high  dynamic  range  without  loss  of  axial  or  lateral 
resolution.  This  was  confirmed  by  estimating  the  correla¬ 
tion  cell  size  based  on  [17]  for  the  quadratic,  images  and 
comparing  with  those  for  B-mode  and  Harmonic,  images 
(Tables  I  and  II).  The  improvement  in  contrast  was  con¬ 
firmed  by  quantitative  measures,  both  on  the  RF  data  and 
the  log-compressed  image  data  after  scan  conversion.  Fur¬ 
thermore,  the  quadratic  images  were  also  shown  to  pre¬ 
serve  the  low-scattering  targets  while  improving  the  UCA 
to  tissue  contrast.  This  was  demonstrated  by  the  computed 
CR  values  (Tables  I  and  II)  and  the  Az  values  from  the 
ROC  curves  shown  in  Fig.  12. 
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ABSTRACT 

It  has  long  been  recognized  that  thermal  lesions  formed  us¬ 
ing  high-intensity  focused  ultrasound  (HIFU)  exhibit  non¬ 
linear  behavior  that  can  be  detected  in  pulse-echo  ultrasound. 
Second  harmonic  imaging  of  freshly  formed  thermal  lesions 
have  consistently  shown  significant  enhancement  in  their  vi¬ 
sualization  confirming  this  nonlinear  behavior.  In  this  pa¬ 
per,  we  describe  a  post-beamforming  nonlinear  filtering  al¬ 
gorithm  based  on  second-order  Volterra  filter  (SVF)  model 
that  separates  the  linear  and  quadratic  components  of  the 
echo  signal  leading  to  significant  enhancement  of  lesion  vi¬ 
sualization.  Images  from  ex  vivo  tissue  samples  are  shown 
to  demonstrate  the  level  of  contrast  enhancement  achieved 
with  the  SVF-based  quadratic  filter  compared  with  standard 
echo  and  2nd  harmonic  imaging  results. 

1.  INTRODUCTION 

Since  early  2001 ,  we  have  used  a  dual-mode  array  described 
in  [1]  to  form  HIFU-induced  thermal  lesions  in  freshly  ex¬ 
cised  degassed  tissue  under  a  variety  of  normal  exposure 
and  over  exposure  conditions.  Single-transmit  focus  im¬ 
ages  were  collected  for  over  100  lesions  before  and  after 
lesion  formation.  These  images  have  consistently  shown  5 
-  7  dB  enhancement  in  the  echogenicity  from  the  lesion  lo¬ 
cation  in  the  standard  echographic  images.  These  results 
were  much  more  consistent  than  the  reported  "flashes"  on 
the  B-scan  images  when  diagnostic  ultrasound  systems  are 
used  to  monitor  HIFU  lesion  formation.  Motivated  by  the 
excellent  investigation  by  P.  P.  Lele  reported  in  [2],  we  hy¬ 
pothesized  that  this  change  in  echogenicity  is  due  to  sta¬ 
ble  microbubbles  that  can  occur  even  at  low  insonaiion  lev¬ 
els  .  Lele  found  that  subhannonic  emission  due  to  microbub- 
bles  showed  a  monotonic  increase  with  intensity  from  150 
mW/cm2  to  1500  W/cm2  without  a  distinct  threshold  for 
emission  (measurements  done  in  vitm  and  in  vivo  at  2.7  and 
1 .8  MHz).  The  consistency  of  the  increase  in  echogenicity 
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at  the  lesion  may  be  explained  by  the  fact  that  the  microbub¬ 
bles  may  already  be  resonant  at  the  imaging  frequency  (same 
as  the  therapeutic  HIFU  beam  when  the  dual-mode  array  is 
is  used),  perhaps  a  result  of  rectified  diffusion. 

The  standard  echographic  images  at  the  fundamental, 
however,  offer  limited  contrast  enhancement  due  to  the  speckle 
phenomenon.  Therefore,  they  could  not  provide  a  reliable 
method  for  mapping  the  boundaries  of  HIFU-induced  le¬ 
sions.  This  lead  us  to  try  to  exploit  the  nonlinear  nature 
of  the  microbubbles  to  enhance  the  visualization  and  map¬ 
ping  of  thermal  lesions.  The  idea  was  that,  if  microbubbles 
are  indeed  present  at  the  lesion  location,  they  will  generate 
nonlinear  echoes  that  may  be  better  suited  for  mapping.  We 
have  initially  investigated  second  harmonic  (SH)  imaging 
as  a  means  of  enhancing  the  lesion  contrast  for  improving 
the  visualization  of  these  images.  SH  images  of  thermal  le¬ 
sions  have  shown  increase  in  the  contrast  on  the  order  of 
22  -  25  dB,  but  with  decreased  dynamic  range  of  the  re¬ 
sulting  images  [1].  A  post-beamforming  nonlinear  com¬ 
pounding  algorithm  was  shown  to  improve  the  contrast  in 
lesion  echogenicity  to  30  -  35  dB  without  loss  in  dynamic 
range  [3].  This  was  achieved  by  compounding  the  funda¬ 
mental  and  the  SH  images  using  spatial  compounding  func¬ 
tions  based  on  the  receive  beamforming  characteristics  of 
the  dual-mode  array  at  the  fundamental  and  the  SH  frequen¬ 
cies.  In  this  paper,  we  describe  a  post-beamforming  filter¬ 
ing  algorithm  for  separating  the  linear  and  quadratic  com¬ 
ponents  of  the  beamformed  data  based  on  the  SVF.  This 
approach  is  superior  to  all  of  the  algorithms  based  on  SH 
imaging.  It  greatly  enhances  lesion  to  tissue  contrast  with¬ 
out  any  loss  in  spatial  resolution.  In  addition,  it  enhances 
the  dynamic  range  of  the  image  thus  greatly  improving  both 
detecting  and  mapping  of  HIFU-induced  lesions,  even  for 
volumetric  lesions. 

2.  IMAGE  FORMATION 

Figure  1  summarizes  the  image  acquisition  and  image  for¬ 
mation  model.  A  64-element  array  optimized  for  maximum 
energy  delivery  at  1  MHz  (Imasonic,  Besan$on,  France) 
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Fig.  1.  Imaging  Setup 


is  used  for  lesion  formation  in  sample  tissue.  Lesions  are 
formed  by  focusing  the  array  at  a  point  within  the  target 
and  maintaining  high-power  output  for  time  intervals  on  the 
order  of  seconds  (1-5  seconds  typical).  The  power  is  in¬ 
terrupted  for  short  intervals  (milliseconds)  to  acquire  image 
data  by  transmitting  short  (pis)  pulse  from  all  64-elements 
and  receiving  on  selected  elements  using  a  matrix  switch. 
Once  the  image  data  set  is  collected,  RF  beamforming  is 
performed  to  form  standard  echographic  images  of  the  tar¬ 
get  region.  Figure  I  also  shows  a  general  post-beamforming 
filter  bank  for  analysis  of  beamformed  data.  For  exam¬ 
ple,  the  filter  bank  could  be  designed  to  perform  harmonic 
decomposition  of  beamfromed  data  (e.g.,  fundamental  and 
2nd  harmonic  imaging).  Alternatively,  the  filter  bank  can 
be  designed  to  separate  the  linear  and  quadratic  (nonlin¬ 
ear)  components  based  on  the  2nd  order  Volterra  model  de¬ 
scribed  below. 


3.  SECOND-ORDER  VOLTERRA  MODEL 

Results  from  [4]  have  shown  the  validity  of  a  second-order 
Volterra  filer  as  a  model  for  pulse-echo  ultrasound  imaging 
data  from  tissue  mimicking  media.  In  this  section,  the  de¬ 
composition  of  received  echo,  i.e.,  output  sequences  only, 
into  linear  and  quadratic  components  by  using  least-squares 
approach  of  second-order  Volterra  model  will  be  considered 
and  the  detail  of  algorithm  implementation  to  pulse-echo  ul¬ 
trasound  imaging  will  be  stated. 


3.1.  Signal  Separation  Model 

The  algorithm  described  in  this  section  is  adapted  from  [5J. 
The  response  of  a  quadratically  nonlinear  system,  y(n  + 1), 
can  be  predicted  by  a  second-order  Volterra  model  of  past 


m  values  as  follows: 

y{n  +  1)  =  yL{n  +  1)  +  yQ{n  +  1) 

m— 1 

y(n  -  i)hL(i) 

i—O 

m-1  m-1 

+  3)yin  *  k)hQ(h  k)  +  e(n),  (1) 

3=0  k=j 

where  hi  (i)  is  linear  filter  coefficients,  hq(j,  k )  represents 
quadratic  filter  coefficients  and  e(n)  is  a  modelling  error 
and/ora  measurement  noise  which  is  assumed  to  be  an  inde¬ 
pendent,  identically  distributed(i.i.d)  random  variable  with 
zero  mean.  That  is,  if  the  model  coefficients  are  known,  the 
echo  signal  can  be  decomposed  into  linear  and  quadratic 
components.  The  latter  can  be  expected  to  better  represent 
the  quadratic  response  of  the  system  than,  say,  the  second 
harmonic  component.  The  model  coefficients  can  be  ob¬ 
tained  by  setting  the  linear  and  quadratic  prediction  problem 
in  Equation  1  to  form  a  set  of  linear  equations.  Recognizing 
that  the  output  is  linear  in  terms  of  the  (unknown)  model 
coefficients,  one  obtains  a  matrix  equation  of  the  form: 

f  =  Gh  +  e,  (2) 

where  the  vector  f,  the  matrix  G  and  the  error  vector  e  are 

f  =  [ y(n  +  l),y(n-h2),...,y(n  +  L)]T 
G=  [y(n),y(n  +  l),...,y(n  +  L-  1)]T 
e  =  [e(n),e(n  4- 1),  ■ . .,e(n  +  L-  1)]T. 

where  the  data  vector,  y,  is  given  by: 

y(n)  =  [y{n),y{n  -  1), y{n  -  2), . . . ,  y(n  -  m  -I- 1), 
V2{n),y{n)y(n  -  1), . . .  ,y2(n  -  m  +  l)]r 
and  the  filter  coefficient  vector,  h,  is  given  by: 

h=  [fri(0),/iz.(l),l»r,(2),...,ht(Tn-  1), 

hQ(0,0),hQ(0,l),...,hQ(m  -  l,m-l)]T. 

The  details  of  the  solution  for  the  coefficients  of  the  SVF 
model  can  be  found  in  [6].  Briefly,  a  minimum-norm  least- 
squares  solution  of  (2)  is  obtained  using  truncated  singular 
value  decomposition,  TSVD.  To  assess  the  performance  of 
the  signal  separation  mode!  in  enhancing  the  lesion  visual¬ 
ization,  we  compute  the  contrast-to-tissue  ratio: 

CTS=10!os“(lfi)  ,3> 

where  ||ygc  lb  and  ||ygr  II2  are  the  l2  norms  of  the  quadratic 
components  from  the  lesion  and  normal  tissue  regions,  re¬ 
spectively.  These  regions  are  easily  identified  under  vari¬ 
ous  imaging  conditions.  For  instance,  for  the  application 
described  in  this  paper,  the  contrast  region  is  the  expected 
location  of  the  thermal  lesion  (often  visible  on  the  standard 
echographic  image). 
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3.2.  Implementation 

The  coefficients  of  the  linear  and  quadratic  components  of 
the  SVF  model  are  obtained  from  the  beamformed  RF  data. 
The  implementation  steps  are  as  follows: 

1 .  Using  the  standard  echographic  image,  select  a  beam- 
formed  RF  data  segment  from  the  expected  lesion  lo¬ 
cation. 

2.  Form  the  linear  systems  of  equations  according  to  (2). 

3.  Define  contrast  region  (within  the  lesion)  and  normal 
tissue  region  for  the  computation  of  the  mean-square 
error  MSE  and  CTR. 

4.  Solve  systems  of  linear  equations  by  using  TS  VD  reg¬ 
ularization  method. 

5.  Apply  second-order  Volterra  filter  to  the  beamformed 
RF  data  throughout  the  pulse-echo  ultrasound  image, 

6.  The  quadratic  component  from  the  SVF  can  be  dis¬ 
played  as  a  separate  image  or  appropriately  compounded 
with  the  linear  component. 

4.  RESULTS  AND  DISCUSSION 

The  experimental  setup  shown  in  Figure  1  was  used  in  ob¬ 
taining  images  of  ex  vivo  tissue  samples  before  and  after 
thermal  lesion  formation  with  a  dual-mode  64-element  ul¬ 
trasound  phased  array  operating  at  1  MHz.  The  standard 
echographic  images  were  formed  by  single  transmit  beams 
[1]  along  the  main  axis  of  the  array  and  dynamic  receive 
focusing  at  each  pixel  in  the  image.  Results  shown  in  [1] 
confirmed  that  echoes  from  the  lesion  location  exhibited  in¬ 
creased  levels  of  second  harmonic  generation  even  at  nor¬ 
mal  exposure  conditions.  Figure  2  shows  the  RF  data  along 
central  axis  of  the  dual-mode  array  before  and  after  lesion 
formation  at  normal  exposure  level  (850  W/cm2  for  4  sec¬ 
onds).  The  figure  also  shows  spectrograms  of  the  RF  data 
(showing  the  frequency  content  of  the  RF  echoes  in  the  ax¬ 
ial  direction).  The  results  show  that  the  echoes  before  le¬ 
sion  formation  are  centered  at  1  MHz  with  no  evidence  of 
2nd  harmonic  component  in  the  tissue  region  (from  80  to 
1 20  mm).  On  the  other  hand,  the  spectrogram  of  the  RF 
data  from  the  same  direction  after  lesion  formation  shows  a 
strong  2nd  harmonic  component  at  the  lesion  location  (90 
-  100  mm).  Figure  3  shows  the  RF  data  along  central  axis 
of  the  dual-mode  array  before  and  after  lesion  formation  at 
normal  exposure  level  (2150  W/cm2  for  3  seconds).  The 
spectrograms  show  that  the  echoes  before  lesion  formation 
are  centered  at  1  MHz  with  small  2nd  harmonic  compo¬ 
nent  in  the  tissue  region  (at  100  mm).  The  spectrogram  of 
the  RF  data  from  the  same  direction  after  lesion  formation 
shows  a  strong  2nd  harmonic  component  starting  at  82  mm. 
This  is  consistent  with  the  shape  of  the  lesion  determined  by 
histological  evaluation  [1]  (tear-drop  cross  section  with  the 
tip  near  the  geometric  center  and  the  base  near  the  front  of 
the  sample).  It  is  interesting  to  note  that,  in  addition  to  the 


increased  2nd  harmonic  generation,  the  echoes  from  the  le¬ 
sion  appear  to  have  wider  bandwidth  compared  to  echoes 
from  the  same  location  before  lesion  formation.  This  is 
reminiscent  of  the  behavior  of  ultrasound  contrast  agents 
(UCAs).  It  provides  indirect  evidence  that  bubble  oscilla¬ 
tion  occurs  during  and  after  lesion  formation  and  lingers  for 
tens  of  seconds  (in  vitro  depending  on  lesion  size). 

In  [1]  we  have  shown  images  of  single  shot  HIFU  le¬ 
sions  at  the  fundamental  and  SH  frequencies  of  the  dual¬ 
mode  transducer.  Both  imaging  modes  produced  satisfac¬ 
tory  images  of  the  lesion  in  these  cases,  with  the  SH  im¬ 
ages  providing  higher  contrast  of  lesion  echogenicity.  How¬ 
ever,  both  of  these  methods  suffer  when  imaging  extended 
lesions.  One  such  example  is  shown  in  Figure  4  before  (left) 
and  after  lesion  formation  at  the  fundamental  (top  pair),  SH 
(middle  pair)  and  using  the  quadratic  component  from  the 
SVF  (bottom  pair),  One  can  easily  observe  the  increase  in 
echogenicity  at  the  lesion  location  (90  -  100  mm  along  the 
axis  of  the  array)  for  all  three  modes.  However,  one  can  see 
that  the  dynamic  range  of  the  fundamental  image  is  limited 
by  the  speckle  pattern  from  tissue  surrounding  the  lesion 
while  the  SH  image  is  limited  by  the  beamforming  artifacts. 
The  lesion  echogenicity  in  this  case  is  nearly  17  dB  above 
the  normal  tissue  echogenicity  in  standard  echographic  im- 
age  (typical  contrast  is  more  like  5  -  7  dB).  The  SH  images 
improve  the  lesion  contrast  to  approximately  25  dB  com¬ 
pared  to  normal  tissue  (22  -  25  dB  typical),  but  beamform¬ 
ing  artifacts  compromise  the  definition  of  the  lesion  bound¬ 
aries  in  the  lateral  direction. 

Using  a  beamformed  echo  data  segment  along  the  main 
axis  of  the  airay,  the  coefficients  of  linear  and  quadratic 
components  of  the  SVF  were  estimated  as  described  in  Sec¬ 
tion  3.2.  The  quadratic  components  of  the  beamformed  RF 
data  were  filtered  out  at  all  pixel  locations.  The  images  be¬ 
fore  and  after  lesion  formation  are  shown  in  Figure  4  at  40 
dB  dynamic  range.  One  can  see  quite  clearly  that  the  grating 
lobe  components  visible  in  the  2nd  harmonic  images  are  ef¬ 
fectively  eliminated  in  the  quadratic  component  images  ob¬ 
tained  through  the  SVF.  In  addition,  the  speckle  components 
conspicuous  in  the  standard  echo  images  is  greatly  reduced. 
This  level  of  enhancement  is  typical  and  has  been  observed 
consistently  in  over  100  experiments  similar  to  the  one  de¬ 
scribed  in  this  paper.  The  reader  can  appreciate  that  the  le¬ 
sion  boundaries  are  well  defined  in  both  the  axial  and  lateral 
directions.  With  the  significant  increase  in  dynamic  range, 
one  can  see  that  both  detection  and  mapping  of  thermal  le¬ 
sions  is  significantly  facilitated  by  the  use  of  the  quadratic 
filter  based  on  the  SVF  model. 

5.  CONCLUSIONS 

Experimental  results  from  ex  vivo  tissue  samples  provide 
the  strongest  evidence  yet  that  thermal  lesions  exhibit  non¬ 
linear  behavior  as  a  propagation  medium.  Using  a  the  SVF 


2002  IEEE  ULTRASONICS  SYMPOSIUM  —  1437 


Fig.  2.  RF  data  and  spectrograms  before(left)  and  after 
(right)  normal  exposure. 
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Fig.  3.  RF  data  and  spectrograms  before  (left)  and  after 
(right)  over  exposure. 


model,  we  have  separated  the  linear  and  quadratic  compo¬ 
nents  in  the  beamformed  RF  echo  data  and  formed  images 
from  the  quadratic  components.  The  quadratic  component 
images  show  significant  enhancement  in  lesion  visualiza¬ 
tion  due  to: 

1.  They  directly  exploit  the  nonlinear  nature  of  freshly 
formed  thermal  lesion  (possibly  due  to  formation  of 
microbubbles). 

2.  Quadratic  component  combines  both  low  frequency 
(close  to  dc)  and  harmonic  frequency  in  forming  non¬ 
linear  echoes.  This  simultaneously  reduces  speckle 
and  beamforming  artifacts  without  loss  in  spatial  res¬ 
olution. 

3.  The  quadratic  kernel  of  the  SVF  rejects  the  additive 
white  Gaussian  noise  components  which  significantly 
improves  the  SNR  of  the  imaging  system  and  enhances 
the  visualization  of  low  echogenicity  regions  in  the 
image. 
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Fig.  4.  Images  before  (left)  and  after  (right)  formation  of 
lesion.  Top:  standard  echo.  Middle:  second  harmonic.  Bot¬ 
tom:  Quadratic. 
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Abstract  -  We  have  recently  demonstrated  the  imaging  capa¬ 
bilities  of  a  prototype  64-element  1  MHz  concave  array  with  100 
mm  radius  of  curvature.  This  array  was  optimized  for  therapeutic 
applications  using  high-intensity  focused  ultrasound  (HIFU).  We 
have  shown  that  this  dual-mode  ultrasound  array  (DMUA)  has  a 
therapeutic  operating  field  ( ThxOF )  that  extends  by  ±3  cm  and 
±2  cm  around  its  geometric  center  in  the  axial  and  lateral  direc¬ 
tions,  respectively.  We  have  also  shown  that  appropriate  apodiza- 
tion  and  accounting  for  element  directivity  along  with  conventional 
synthetic  aperture  beamforming  produce  a  50  dB  imaging  field  of 
view  ( IxFOV )  larger  than  the  ThxOF.  In  addition,  the  spatial  reg¬ 
istration  of  imaging  targets  is  as  accurate  as  commercially  avail¬ 
able  scanners.  In  this  paper,  we  present  results  from  an  image- 
based  refocusing  algorithm  whereby  images  formed  by  the  DMUA 
are  used  to  identify  a  refocusing  target  and  a  set  of  critical  points 
where  the  incident  power  is  to  be  minimized.  The  algorithm  is  val¬ 
idated  experimentally  in  tissue  mimicking  phantom  with  strongly 
scattering  ribs  placed  between  the  DMUA  and  the  target.  These 
results  demonstrate  what  is  potentially  the  most  powerful  advan¬ 
tage  of  the  use  of  DMUAs  in  image-guided  surgery.  Namely,  the 
inherent  registration  between  the  imaging  and  therapeutic  coordi¬ 
nate  systems.  This  allows  for  direct  definition  of  targets  and  any 
surrounding  critical  structures  to  be  avoided  to  minimize  the  col¬ 
lateral  damage.  With  these  capabilities,  DMUAs  may  provide  a 
most  powerful  paradigm  for  image-guided  surgery. 

I.  INTRODUCTION 

Recently,  non-invasive  and  minimally  invasive  thermal 
surgeries  have  become  widely  accepted  in  clinics  and  hos¬ 
pitals  worldwide.  Reliable  imaging  techniques  for  aiming, 
monitoring  and  visualization  can  help  clinics  gradually  ac¬ 
cept  non-invasive  and  minimally  invasive  thermal  surgeries. 
Image  guidance  methods  based  on  MRI  [1],  CT  [2]  and  ul¬ 
trasound  [3]  have  been  recently  proposed. 

Most  of  the  image  guidance  systems,  including  some  ul¬ 
trasound  systems,  have  separate  image  guidance  sub-system 
and  therapeutic  sub-system.  The  advent  of  piezocomposite 
transducer  technology  has  made  it  feasible  to  design  and  fab¬ 
ricate  dual-mode  arrays  capable  of  high-power  therapeutic 
delivery  while  having  a  wide  enough  bandwidth  for  pulse- 
echo  operation.  One  such  DMUA  has  been  described  and 
introduced  in  our  previous  reports  (e.g.  [3]).  The  ThxOF 
for  this  prototype  was  experimentally  shown  to  match  the 
theoretically  predicted  profile  based  on  array  and  element 
directivity  analysis  [3],  Furthermore,  the  IxFOV  was  also 


shown  to  extend  beyond  the  ThxOF  for  this  prototype  [4]. 

The  main  advantage  of  this  approach  is  the  inherent  regis¬ 
tration  between  the  therapeutic  and  imaging  coordinate  sys¬ 
tems.  Once  the  real-time  imaging  capability  is  available  for 
DMUAs,  one  can  envision  a  treatment  paradigm  in  which 
the  HIFU  beam  is  guided  in  real-time  based  on  imaging  of 
the  treatment  field  by  the  same  transducer.  The  doctor  only 
needs  to  locate  the  target  on  the  image  and  the  system  auto¬ 
matically  generates  the  amplitude  and  phase  distribution  to 
focus  the  HIFU  pulse  at  the  desired  target.  In  addition,  it 
is  possible  to  define  additional  points  in  the  treatment  field 
where  the  incident  power  is  to  be  minimized  to  reduce  or 
eliminate  collateral  damage,  e.g.  to  nearby  bone  structures. 
In  this  paper,  we  present  the  first  experimental  verification 
of  an  image-based  refocusing  algorithm  employing  images 
from  the  DMUA  prototype  to  identify  the  refocusing  target 
as  well  as  a  number  of  critical  structures  to  be  avoided.  One 
obvious  application  of  this  capability  is  targeting  tumors  in 
the  liver  that  are  partially  obstructed  by  the  ribcage  [5],  This 
capability  is  unique  to  DMUAs  due  to  the  inherent  regis¬ 
tration  between  the  imaging  and  therapeutic  coordinate  sys¬ 
tems. 

II.  Theory 

A.  Synthetic  Aperture  (SA)  Imaging 

Images  were  obtained  using  a  full  synthetic  aperture  tech¬ 
nique  [6].  The  image  pixel  at  coordinates  (xp,zp)  was  there¬ 
fore  computed  by  [3]: 

64  64 

I{xp,Zp)  —  ^  ^  'y  ^  Aj  •  13 j  •  SiJ[(Rip  A  Rjp)/c\,  (1) 

i=l  j= 1 

where  i  is  the  transmit  element  index,  j  is  the  receive  ele¬ 
ment  index.  Af  is  the  transmit  apodization  weight  at  element 
i,  Bj  is  the  receive  apodization  weight  at  element  j,  R-tp  and 
RjP  are  the  distance  from  the  transmit  and  receive  element, 
respectively,  from  the  image  pixel,  c  is  the  speed  of  sound  in 
the  medium  being  imaged,  and  Si,j(t)  is  the  echo  acquired 
where  transmitting  with  element  i  and  receiving  with  ele¬ 
ment  j.  Specialized  image  reconstruction  programs  were 
written  to  evaluate  the  intensity  at  ( xp ,  zp). 

B.  Single-Transmit  Focus  (STF)  Imaging 

The  synthetic  aperture  imaging  algorithm  described  above 
can  be  used  to  produce  higher  quality  conventional  images 
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than  can  be  expected  from  a  given  array.  However,  due  to  the 
fact  that  transmit  beams  are  synthesized  by  superposition  of 
single-element  transmit  patterns,  the  nonlinear  interactions 
of  the  real-time  transmit  beams  cannot  be  accounted  for  us¬ 
ing  the  imaging  algorithms.  Therefore,  we  have  modified 
our  64-channel  phased  array  driver  to  allow  for  pulsed  trans¬ 
mission  on  all  64  channels  simultaneously.  This  allowed  us 
to  use  the  full  power  of  the  transmit  beams  and,  therefore, 
observe  their  nonlinear  interactions  with  the  tissue  media. 
The  image  formation  process  due  to  a  single  transmit  fo¬ 
cused  beam  is  a  modified  version  of  Eq.  (1)  as  follows: 

64 

I(xp,  zp)  =  ]T  sA(Ro  +  Rir)/C 1’  (2) 

j= i 

where  s j  ( t )  is  the  received  waveform  at  element  j  due  to  the 
transmitted  beam  and  Rq  is  a  fixed  distance  determined  by 
the  focal  depth  of  the  transmit  beam.  All  other  quantities  in 
Eq.(2)  are  the  same  as  their  counterparts  in  Eq.  (1). 


C.  Special  Considerations 

In  conventional  phased-array  beamforming,  grating-lobe 
levels  are  determined  by  the  array  pattern  as  well  as  the 
element  directivity.  Both  of  these  factors  are  important  in 
DMUAs  as  the  array  elements  typically  have  high  directiv¬ 
ity. 

1)  Array  Gain  Compensation  (GC)  Algorithm:  The  in¬ 
tensity  gain  is  defined  as  the  peak  field  intensity  in  W /cm2 
divided  by  the  average  W / cm2  output  from  the  array  at  a  set 
of  field  locations.  In  order  to  correctly  calculate  the  contrast 
ratio  of  the  image  in  a  speckle  region,  gain  compensation 
is  required.  Since  the  array  is  dynamically  focused  on  both 
transmit  and  receive,  2-way  compensation  is  utilized. 

2)  Element  Directivity  (ED)  Algorithm:  Accounting  for 
element  directivity  in  the  beamforming  algorithm  reduces 
the  grating  lobes  and,  consequently,  improves  the  dynamic 
range.  Furthermore,  accounting  for  element  directivity  af¬ 
fects  the  SNR  of  the  resulting  image.  The  element  directiv¬ 
ity  weight  function  is: 


D(6pe)  = 


sin[fcw  sin(0pe)/2] 
kw  sin(0pe)/2 


(3) 


where  k  —  A  is  the  wavelength,  w  is  the  element  width, 
9pe  is  the  angle  between  the  axis  of  the  radiating  (receiving) 
element  and  the  vector  from  the  element  center  to  the  pixel. 


D.  Refocusing  in  the  Presence  of  Obstacles 
Assume  we  have  an  Ar-element  array  with  arbitrary,  but 
known,  geometry  radiating  into  a  homogeneous  half  space. 
The  objective  is  to  maximize  the  array  intensity  gain  at  a 
target  point  fix  while  minimizing  the  incident  power  at  a  set 
of  critical  points,  fic(t),  *  —  1,2,...,  Me-  This  is  an  opti¬ 
mization  problem  that  can  be  solved  using  Lagrange  multi¬ 
pliers  or  a  regularized  minimum-norm  least  squares  solution 
as  follows  [5]: 


•  Let  the  raw  vector  h  =  [hi(fT),h-2{ff), . . . , /ijv'(fix)] 
be  the  array  directivity  vector  at  fix  (note  that  /u-(Lr) 
is  the  directivity  of  the  fcth  element  at  fix). 

•  Form  the  matrix  He  from  the  raw  vectors  h,  = 
[hi  (fc(i)),h2(fc(i)),...,  hN  (fib  (*))]. 

•  Form  the  weighting  matrix,  Wc  =  (HcH£  +  7I)-1 

where  I  is  the  identity  matrix  and  7  is  an  appropriately 
chosen  regularization  parameter. 

•  The  optimal  complex  array  driving  vector  is  given  by 
u  =  Wch£  (hxWch^)'1. 

Both  the  target  and  the  critical  points  can  be  derived  from 
images  formed  by  the  DMUA  in  synthetic  aperture  or  in 
single-transmit  focus  modes.  This  is  probably  the  most  pow¬ 
erful  feature  of  the  use  of  the  DMUA  in  image-guided  appli¬ 
cations. 

III.  Experiment  Setup 

Figure  1  (left)  shows  the  image  acquisition  model.  A  64- 
element  concave  linear  phased  array  is  used  for  the  image 
acquisition.  The  array  can  be  connected  to  a  64-channel 
driver  that  can  operate  in  CW  therapeutic  mode  or  pulsed 
imaging  mode.  It  can  also  be  connected  to  a  standard  pulser 
receiver  through  a  matrix  switch.  Synthetic  aperture  images 
were  obtained  with  the  latter  approach  while  single  trans¬ 
mit  images  were  obtained  with  the  64-channel  transmitter. 
Figure  1  (right)  shows  a  schematic  of  the  DMUA  with  an 
imaging  target  containing  three  strongly  scattering  ribs  on 
top  of  a  tissue-mimicking  phantom  containing  a  strong  scat- 
terer  near  the  geometric  center  of  the  DMUA.  The  target  is 
a  thermocouple  wire  0.1  mm  in  diameter  and  the  ribs  are 
plastic  bars  with  diameter  of  8  mm.  They  are  positioned  in 
a  plane  just  over  40  mm  from  the  vertex  of  the  DMUA  with 
approximate  spacing  of  20  mm  between  them. 

A  modified  Technos  MPX  system  from  ESAOTE,  Genoa, 
Italy  is  used  to  acquire  images  for  comparison  with  images 
of  the  same  target  from  the  DMUA.  The  system  is  modi¬ 
fied  to  allow  the  upload  of  high-quality  beamformed  RF  data 
from  the  Technos  system  to  computer.  A  CA  421  convex 
abdominal  probe  is  used  to  acquire  images  from  the  same 
target  as  the  DMUA. 

IV.  Results 

A.  Synthetic  Aperture  Imaging 

Two  examples  of  SA  imaging  with  the  DMUA  are  shown 
to  demonstrate  its  imaging  capability.  The  first  is  the  result 
of  a  3D  contrast  target  in  a  quality  assurance  phantom  and 
the  second  is  a  tissue  mimicking  phantom  with  rib  obstacles 
used  in  the  refocusing  experiments  described  below. 

1)  3D  Contrast  Phantom:  Figure  2  shows  grayscale  im¬ 
ages  (50  dB)  of  an  egg-shaped  3D  contrast  object  within 
the  CIRS  Model  55  3D  contrast  phantom.  The  image  on 
the  left  was  obtained  using  the  DMUA  while  the  one  on  the 
right  was  obtained  using  the  Technos  with  the  CA421  con¬ 
vex  probe.  The  contrast  ratio  (CR)  for  the  egg-shaped  con- 
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Figure  1 :  Dual-mode  system  for  lesion  and  image  formation  and 
imaging  target  with  obstacles  used  in  refocusing  experiments. 
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Figure  2:  Images  (50  dB)  from  the  CIRS  Model  55  3D  contrast 
phantom.  Left:  DMUA  w/  GC;  Right:  Technos. 

trast  target  was  computed  as  follows: 

where  Ic[i)  [/*(!)]  is  the  average  intensity  in  a  3  x  3  mm2  re¬ 
gion  inside  [outside]  the  contrast  target  and  Nr  is  the  number 
of  regions.  The  CR  value  determined  from  image  of  DMUA 
(with  gain  compensation)  and  the  image  from  the  CA  421 
probe  are  8.35  dB  and  9.12  dB,  respectively.  These  values 
are  quite  consistent  with  the  manufacturer  specified  contrast 
of  9  ±  0.5  dB.  It  should  be  mentioned  that,  without  com¬ 
pensating  for  intensity  gain  and  element  directivity  for  the 
DMUA,  the  contrast  ratio  would  have  been  3.3  dB  instead 
of  8.35  dB.  It  is  also  interesting  to  note  that  the  distal  edge 
of  the  egg-shaped  object  marks  the  far  end  of  the  50  dB  Ix- 
FOV  of  the  DMUA.  One  can  visually  observe  the  loss  of 
contrast  in  that  region.  Spatial  and  contrast  resolution  re¬ 
sults  for  various  beamforming  conditions  are  shown  in  [4]. 

2)  Tissue-Mimicking  Phantom  with  Obstacles:  Figure  3 
shows  a  grayscale  image  (50  dB)  of  the  phantom  shown  in 
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Figure  3:  Grayscale  image  of  the  tissue-mimicking  phantom  with 
ribs  shown  in  Figure  It  acquired  using  the  DMUA  in  synthetic 
aperture  mode  (50  dB). 

Figure  1 .  The  image  was  obtained  using  the  full  synthetic 
aperture  technique  described  in  Section  II-A.  The  result  is 
consistent  with  the  schematic  with  good  spatial  resolution 
as  well  as  contrast.  In  addition  to  the  target  and  the  ribs,  the 
front  and  back  of  the  tissue-mimicking  phantom  can  be  eas¬ 
ily  discerned  (at  80  and  1 20  mm,  respectively).  Even  though 
the  ribs  may  be  outside  the  50  dB  IxFOV  of  the  DMUA,  they 
are  sufficiently  strong  to  standout  giving  strong  specular  re¬ 
flections.  Furthermore,  the  ribs  are  all  spatially  registered 
correctly  without  distortion. 

B.  Single-Transmit  Focus  Imaging 

Single-transmit  focus  imaging  can  be  extremely  useful 
in  conjunction  with  DMUAs  in  image-guided  surgery.  By 
imaging  the  target  tissue  using  the  delay  and  magnitude  pro¬ 
file  of  the  therapeutic  beam,  it  is  possible  to  visualize  any 
strongly  scattering  structures  directly  affected  by  this  beam. 
Figure  4  shows  the  STF  image  obtained  using  a  transmit  fo¬ 
cus  at  0  lateral  and  100  mm  axial  (geometric  focus).  The  im¬ 
age  is  displayed  with  40  dB  dynamic  range  and  shows  that 
two  of  the  three  ribs  are  visible  in  addition  to  the  target  near 
the  geometric  center.  This  result  implies  that  the  transmit 
focus,  when  used  in  therapeutic  mode,  may  cause  excessive 
heating  at  one  or  two  of  the  ribs  visible  on  the  image. 

While  the  STF  image  clearly  has  lower  contrast  than  the 
full  SA  image  shown  in  Figure  3,  it  still  shows  the  target  and 
two  of  the  ribs  with  good  spatial  registration.  Therefore,  it  is 
still  useful  for  identifying  targets  and  potential  obstacles  as 
described  in  the  next  subsection. 

C.  Refocusing  in  the  Presence  of  Obstacles 

The  image  shown  in  Figure  4  was  used  to  identify  the  co¬ 
ordinates  of  the  visible  ribs  as  well  as  define  a  target,  point 
for  refocusing  within  the  tissue  mimicking  phantom.  The  40 
dB  grayscale  images  shown  in  Figure  5  show  the  results  of 
refocusing  without  (left)  and  with  (right)  taking  the  rib  loca¬ 
tion  into  consideration.  Both  of  these  images  where  obtained 
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Figure  4:  STF  image  (40  dB)  of  the  target  shown  in  Figure  3. 
Strong  echoes  from  the  three  ribs  are  evident  while  the  needle  is 
less  visible. 
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Figure  5:  STF  images  (40  dB)  using  refocusing  without  (left)  and 
with  (right)  taking  the  critical  points  are  the  ribs  into  consideration. 

with  a  transmit  beam  focused  at  the  defined  target  location 
at  5  mm  lateral  and  100  mm  axial.  However,  the  image  on 
the  left  was  obtained  using  the  algorithm  described  in  Sec¬ 
tion  II-D.  One  can  clearly  see  the  increased  visibility  of  the 
target  compared  to  the  visibility  of  the  ribs.  This  is  due  to 
the  reduction  of  incident  power  at  the  ribs  as  predicted  by 
the  weighting  algorithm  described  in  Section  II-D. 

V.  Discussion  and  Conclusions 

Characterization  of  the  image  quality  from  the  DMUA 
demonstrates  the  feasibility  of  using  this  array  in  image- 
guided  surgical  applications.  The  spatial  and  contrast  res¬ 
olutions  of  this  array,  while  not  sufficient  for  high-quality 
diagnostic  applications,  are  consistent  with  the  bandwidth 
and  geometry  of  the  array.  Furthermore,  the  speckle  com¬ 
ponents  of  the  beamformed  data  is  suitable  for  relevant  sig¬ 
nal  processing  algorithms  such  as  displacement  tracking  and 
temperature  estimation. 

This  paper  has  also  presented  the  first  experimental  verifi¬ 
cation  of  the  feasibility  of  using  image-based  refocusing  for 
targeting  tissue  structures  while  avoiding  critical  obstacles 


such  as  the  ribs.  This  may  be  important  if  HIFU  arrays  are 
to  be  used  for  targeting,  for  example,  liver  tumors  that  may 
be  partially  obstructed  by  the  rib  cage.  This  result  is  not  to  be 
understood  as  a  form  of  aberration  correction  since  the  refo¬ 
cusing  algorithm  assumes  knowledge  of  the  array  directivity 
at  the  refocusing  target  and  the  critical  points.  Nonetheless, 
since  these  same  assumptions  are  used  in  beamforming  the 
image  in  the  first  place,  the  algorithm  is  quite  robust  against 
distortions  due  to  inaccuracies  in  the  speed  of  sound  and 
other  tissue  properties.  The  only  requirement  is  that  the  SA 
or  the  STF  image  used  in  refocusing  provides  a  recognizable 
map  of  the  treatment  region. 

The  image  quality  currently  available  for  DMUAs  with 
image  formation  techniques  described  herein  may  not  be 
sufficient  to  justify  identifying  more  complex  lower  scatter¬ 
ing  critical  structures.  Furthermore,  we  are  not  suggesting 
that  targets  such  as  tumors  can  be  recognized  on  images  pro¬ 
vided  using  the  DMUA.  However,  continued  improvement 
of  piezocomposite  transducer  technologies,  especially  im¬ 
proved  bandwidths  of  high-power  arrays,  may  bring  about 
sufficient  improvement  in  the  image  quality  to  make  that 
possible.  In  any  event,  the  current  image  quality  may  allow 
for  recognizing  landmarks  on  the  image  that,  with  the  help 
of  more  refined  imaging  systems,  can  be  used  in  the  target¬ 
ing  of  tumors  and  less  recognizable  tissue  structures  as  they 
appear  on  DMUA  images. 
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ABSTRACT 

In  this  paper,  we  demonstrate  that  a  second-order  Volterra 
filter  (SVF)  can  adequately  model  pulse-echo  imaging  data 
acquired  from  tissue  media.  Furthermore,  we  demonstrate 
that  the  coefficients  of  both  the  linear  impulse  response  and 
the  second-order  kernel  can  be  robustly  estimated  from  imag¬ 
ing  data  using  pseudorandom  binary  input  sequences.  The 
validation  of  the  SVF  as  an  appropriate  model  for  ultra¬ 
sonic  imaging  will  allow  for  more  intelligent  use  of  post¬ 
beamforming  signal  processing  for  enhancing  image  con¬ 
trast  and  producing  quantitative  imaging. 

1.  INTRODUCTION 

Modem  ultrasonic  imaging  is  moving  rapidly  towards  new 
paradigms  employing  post-beamforming  filtering  combined 
with  nonlinear  imaging  modes  [9, 1],  Native  tissue-harmonic 
(NH)  and  contrast-agent  (CA)  imaging  are  currently  the  lead¬ 
ing  applications  being  investigated  in  many  laboratories  around 
the  world  [4, 1, 5].  Figure  1  represents  a  general  scheme  of  a 
modem  ultrasound  data  acquisition  and  reconstruction  em¬ 
ploying  post-beamforming  filter  bank  for  signal  separation 
and  compounding.  Once  the  image  data  set  is  collected,  RF 
beamforming  is  performed  to  form  standard  echographic 
images  of  the  target  region.  Post-beamforming  filterbank  is 
used  to  separate  signal  components  of  beamformed  RF  data 
based  on  space- time  processing  of  echoes  from  different  di¬ 
rections  [8]  or  appropriate  nonlinear  filtering  depending  on 
a  known  or  identified  model  of  the  nonlinear  system.  As 
an  example  of  the  latter,  nonlinear  imaging  used  in  prac¬ 
tice  on  medical  scanners  send  transmit  pulses  centered  at  a 
frequency  /0  and  detects  signal  components  at  2/0  to  pro¬ 
duced  images  with  improved  contrast  with  or  without  con¬ 
trast  agents.  Until  very  recently,  this  was  the  approach  taken 
for  both  NH  and  CA  imaging.  A  better  understanding  of  the 
nature  of  nonlinearity  as  well  as  the  availability  of  advanced 
VLSI  for  signal  processing  have  been  the  driving  forces  be¬ 
hind  the  development  of  more  sophisticated  approaches  to 

Funded  in  part  by  Grant  DAMD  17-01-1-330  from  the  US  Aimy  Med¬ 
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post-beamforming  signal  processing. 

In  this  paper,  we  demonstrate  the  applicability  of  the  SVF  to 
the  ultrasound  pulse-echo  imaging  problem.  Furthermore, 
we  demonstrate  that  the  SVF  coefficients  can  be  robustly  es¬ 
timated  from  pulse-echo  data,  even  in  the  presence  of  strong 
speckle  noise. 

_ * _ 


Fig.  1.  Image  Acquisition  and  Reconstruction  Algorithm 

2.  IMAGE  FORMATION  MODEL 
Synthetic  Aperture  Imaging  [6]:  The  image  pixel  at  coor¬ 
dinates  (xr)  zp)  is  computedby  (Fig.  2): 

IV  N 

I(xp,Zp)  =  Ai  ■  Bj  •  Sij  [(Rjp  +  Rjp)  /c]  ,  (1) 

i=l  3  =  1 

where  i  is  the  transmit  element  index,  j  is  the  receive  el¬ 
ement  index,  A<  is  the  transmit  apodization  weight  at  ele¬ 
ment  t,  Bj  is  the  receive  apodization  weight  at  element  j, 
Rip  and  Rjp  are  the  distances  from  the  transmit  and  receive 
elements,  respectively,  from  the  image  pixel,  c  is  the  speed 
of  sound  in  the  medium  being  imaged,  and  «y(t)  is  the 
echo  acquired  when  transmitting  with  element  i  and  receiv¬ 
ing  with  element  j. 
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Fig.  2.  Coordinate  system  used  in  the  synthetic  aperture 
imaging  system. 


Post  Beamforming  Filtering:  Since  beamforming  is  per¬ 
formed  in  the  RF  domain,  all  the  frequency  components  in 
the  received  echo  signals  are  retained.  This  is  important 
since  the  echo  signals  can  be  expected  to  contain  a  mix  of 
harmonics  (and  possibly  subharmonics)  depending  on  the 
nonlinearity  of  the  tissue  medium  being  imaged.  This  is 
especially  true  for  the  single-transmit  imaging  mode  where 
the  transmit  power  is  sufficiently  high  to  produce  harmonic 
components  in  nonlinear  tissue  media.  Postbeamforming 
filtering  can  be  used  to  isolate  specific  harmonic  compo¬ 
nents  and/or  enhance  axial  resolution  if  the  SNR  of  the  sys¬ 
tem  is  sufficiently  high.  Algorithms  for  postbeamforming 
filterbank  image  reconstruction  for  pulse-echo  ultrasound 
can  be  found  in  [9]. 

2.1.  Filter  Bank  Design 

2.1.1.  Linear  Imaging  Model  for  Coded-Excitation 

We  begin  with  an  overview  of  the  proposed  coded-excitation 
imaging  modality  [8,  9].  We  assume  that  the  scatterers  in 
the  region  of  interest  (ROI)  are  located  on  a  uniform  sec¬ 
tor  (r,  4>)  grid.  A  discretized  model  for  the  receive  signal, 
denoted  as  f(k),  can  be  written  as 

N*  Nr 

/(*0  =  S  s^i)9iik  ~  i )  +  «(*)»  (2) 

<=t  j=‘ 

where  Nr  is  the  grid  size  determined  by  the  range  interval 
size  and  the  sampling  rate,  n(k)  is  the  noise,  gi{k-j)  is  the 
system  impulse  response  created  by  a  unit  strength  scatterer 
located  on  the  at  the  j-th  range  grid  position  along  the  angle 
<j> i,  and  Si(j)  is  the  actual  scatterer  strength  at  that  location. 

Let  G,  =  [ffi.i  ,9il2,  • .  ■  and  g(J  be  the  vector 

version  of  -  j).  Then  Equation  2  can  be  rewritten  as 

/  =  Gs  +  n,  (3) 

where  /,  a  and  n  are  the  vector  versions  of  f(k),  s{k )  and 
n(k),  respectively.  The  imaging  operator  G,  defined  as 


norm  least  square  estimate  (denoted  by  a  below)  of  the  scat¬ 
terer  vector  a  can  be  expressed  as 

s  =  G,t(GG*t)'f,  '  (5) 

where mt  denotes  the  Hermitian  transpose  of  a  matrix.  This 
is  the  familiar  minimum-norm  least-square  solution  to  the 
image  reconstruction  problem.  We  have  shown  previously 
[9]  that,  under  realistic  assumptions,  a  singular  value  de¬ 
composition  of  the  psudoinverse  can  be  obtained  using  the 
discrete  Fourier  transform  (DFT).  It  is  important  to  note 
here  that,  while  the  solution  given  in  (5)  is  based  on  a  super¬ 
position,  it  can  be  generalized  to  the  nonlinear  case  when 
the  appropriate  mode!  is  used.  The  following  section  de¬ 
scribes  a  commonly  used  nonlinear  wave  propagation  model 
used  for  this  purpose. 

2.1.2.  Nonlinear  Model 

The  Khokhlov-Zabolotskaya-Kuznetsov(KZK)  equation  has 
been  validated  analytically  and  experimentally  as  an  appro¬ 
priate  model  for.finite-amplitude  (nonlinear)  wave  propaga¬ 
tion  in  soft  (water  or  tissue)  media.  It  is  typically  solved 
numerically  in  normalized  form  in  terms  of  the  pressure,  p, 
and  the  coordinates,  (x,  t).  For  example,  for  a  circularly 
symmetric  source  with  radius,  a  and  radius  of  curvature,  d , 
the  following  normalized  variables  are  defined  in  cylindrical 
coordinates: 


P  =  p/po,  a  —  z/d,  p  —  r/a,  r  =  Uot' 


where  po  is  the  hydrostatic  pressure  and  tu0  is  the  center 
frequency  of  the  excitation  signal.  The  normalized  KZK 
equation  takes  the  form: 


8P 

Oc 


J_  fT  / 18P\ 


.  .  .82P  xr  rfiP 

dr+AW+NPW 

(6) 


In  this  equation,  the  following  parameter  are  defined: 


r  _  £o  tup  a2 
d  2cod’ 


A  =  atjd  = 


flw  pd 

~2qjT 


dpwppo 

Po<% 


where  G  is  a  focusing  gain,  A  is  an  absorption  term,  and 
N  is  the  nonlinearity  parameter  of  the  medium.  A  complete 
description  of  the  KZK  is  beyond  the  scope  of  this  paper. 
For  our  purposes,  however,  it  suffices  to  say  that  this  was 
solved  numerically  and  validated  experimentally  in  differ¬ 
ent  media.  We  use  this  equation  to  test  the  validity  of  the 
SVF  model  for  nonlinear  ultrasonic  wave  propagation  and 
the  system  identification  algorithm  used  to  extract  the  SVF 
coefficients. 


G  =  [G1,G2 . G/v*],  (4) 

is  a  (typically,  underdetermined)  matrix  operator  consisting 
of  N$  linear  convolution  matrices.  Thus,  a  unique  minimum- 


2.1.3.  The  Volterra  Filter 

Nonlinear  processes  can  be  modelled  by  general  p-th  order 
Volterra  system  [7],  Due  to  nature  of  the  problem  and  com¬ 
putational  complexity,  we  consider  only  the  second-order 
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Volterra  filter.  The  output  of  a  SVF  can  be  given  by  : 

•  N- 1  N-l  ' 

y(ri)  =  fto+^3  h(k)x(n-k)+  ^  «(*. j>(n-» )x(n-j) 
fc=o  i,j= o 

(7) 

where  h(k)  represents  linear  impulse  response  and  q(m,n) 
is  quadratic  kernel.  The  quadratic  kernel  is  assumed  to  he 
symmetric,  i.e.,  q{i,j)  =  q(j,i)  and  ho  is  assumed  to  be 
zero.  A  parallel  realization  of  the  SVF  is  shown  in  Fig.  3. 
The  identification  problem  in  this  paper  is  to  find  linear  im- 


Fig.  3.  Parallel  realization  of  a  second-order  Volterra  filter. 

pulse  response,  h(n)  and  quadratic  kernel,  q(i,j)  when  in¬ 
put  and  output  are  known.  Furthermore,  we  restrict  the  input 
to  a  class  of  psudorandom  binary  sequences  (e.g.,  Barker 
codes  or  maximal  length  sequences)  that  are  readily  avail¬ 
able  in  many  ultrasound  and  RADAR  imaging  systems. 

2.1.4.  Identification  of  Quadratic  Nonlinear  Systems 

Identification  of  quadratic  nonlinear  system  algorithm  in  this 
section  is  based  on  a  digital  method  of  modelling  quadrati- 
cally  nonlinear  systems  with  general  random  input  discussed 
in  [2].  Equation  7  can  be  represented  in  the  frequency  do¬ 
main  in  term  of  the  DFT  of  y{n): 

M 

Y(fk)  =  H(fk)X(fk)  +  £  QUi,fi)X(fi)X(f,) 

U—  (M-l) 
i+j—k 

(8) 

where  linear  transfer  function,  H (/*),  and  quadratic  trans¬ 
fer  function,  Q(fhfj),  are  unknown.  M  is  the  length  of 
input  sequence  divided  by  two.  Kim  and  Powers  [2]  pro¬ 
posed  a  general  DFT-based  solution  for  the  system  identi¬ 
fication  problem  using  higher-order  spectra.  We  skip  the 
details  of  the  algorithm,  but  state  that  (8)  can  be  written  in 
matrix  form: 

Y{k)  =  x*h,  (9) 

where  the  DFT  index,  k,  denotes  the  frequency  /*  and  the 
vector  h  contains  the  coefficients  H (k)  and  all  Q(i,j)  ■ 
i  +  j  =  k.  To  solve  the  variable  h,  both  sides  of  (9)  are 
multiplied  from  the  left  by  x*  and  the  expectation  operator 
is  applied  resulting  in, 

JS[x*y(fc)]  =  £[**x‘]h.  (10) 


Equation  (10)  is  linear  system  of  equations,  which  can  be 
solved  uniquely  if  {■E[x*x‘]}-1  exists.  The  vector  h  is 
given  by: 

h=  {£[x*x‘]}-1£[x*y(fc)]  (11) 

To  solve  for  every  coefficient  in  the  linear  transfer  function 
and  quadratic  transfer  function,  (1 1)  is  solved  repeatedly  for 
k  =  0, . . . ,  M.  For  quadratic  transfer  function,  the  un¬ 
known  variables  are  solved  along  diagonal  lines  i  +  j  = 
k  in  the  2D  frequency  plane.  In  addition,  the  symmetry 
properties  of  g(m,  n)  are  exploited  to  obtain  a  computa¬ 
tionally  efficient  solution.  Equation  1 1  was  solved  for  the 
case  of  binary  pseudorandom  input  sequences.  The  orig¬ 
inal  algorithm  calls  for  a  full  set  of  M/2  independent  in¬ 
put  sequences  to  insure  a  stable  inverse.  However,  we  have 
tested  the  use  of  a  regularized  solution  to  (11)  using  the 
pseudoinverse  with  a  minimum  set  of  training  sequences. 
Robust  estimates  of  the  SVF  coefficients  were  obtained  us¬ 
ing  a  small  set  of  training  sequences.  This  is  important  for 
a  variety  of  imaging  scenarios.  For  instance,  high  framer- 
ate  imaging  where  we  cannot  afford  to  send  tens  of  training 
sequences.  Another  example  is  imaging  microbubble  con¬ 
trast  agents  where  the  time-invariance  of  the  system  holds 
for  a  few  milliseconds  requiring  system  identification  to  be 
completed  using  a  few  transmisions. 

3.  RESULTS  AND  DISCUSSION 
Simulations:  Figure  4  shows  the  result  of  simulating  the 
propagation  of  a  3-cycle  sinusoidal  pulse  using  the  KZK 
equation  (solid  lines)  in  time  and  frequency  domain.  The 
result  show  a  classical  nonlinear  distortion  of  a  narrowband 
pulse  where  the  negative  pressure  peaks  begin  to  distort  to 
produce  sharper  negative  to  positive  transition.  If  the  pulse 
is  allowed  to  propagate  sufficiently  long,  the  distortion  pro¬ 
duces  a  sawtooth  waveform  and  eventually  leads  to  a  shock 
wave. 

The  SVF  coefficients  were  obtained  using  the  system  iden¬ 
tification  algorithm  described  in  Section  2.1.4  with  64-chip 
maximal  length  sequences.  The  predicted  outputs  from  the 
SVF  (dashed-dotted),  the  linear  impulse  response  (dotted), 
and  the  quadratic  kernel  (dashed)  are  also  shown  for  com¬ 
parison,  It  can  be  seen  that  the  linear  impulse  response 
predicts  the  linear  output  in  the  fundamental  energy  band 
(around  3  MHz).  On  the  other  hand,  the  quadratic  kernel 
captures  energy  primarily  from  the  2nd  harmonic  and  low 
frequency  bands  corresponding  to  second  order  nonlinear¬ 
ity.  It  should  be  noted  that  the  quadratic  kernel  does  not 
capture  the  3rd  harmonic  band  clearly  present  in  the  KZK 
output.  This  is  a  potential  limitation  of  the  SVF  approach 
when  the  system  to  be  identified  is  highly  nonlinear.  The 
simulation  shown  here,  however,  is  an  extreme  case  where 
the  nonlinearity  is  higher  that  can  be  expected  in  practice. 
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Fig.  4.  Propagation  of  a  3-cycle  sinusoidal  pulse  using  the 
KZK  model. 


Experiment:  Figure  5  shows  the  received  signal  from  a 
3-cycle  sinusoidal  pulse  into  a  tissue  mimicking  phantom 
used  for  quality  assurance  testing  of  medical  imaging  scan¬ 
ners.  The  strong  reflection  at  6  cm  is  due  to  a  wire  resolution 
target  positioned  within  a  speckle  region  to  simulate  reflec¬ 
tions  from  human  liver.  We  used  64-chip. maximal  length 
sequences  for  identifying  the  SVF  coefficients  near  the  wire 
target.  For  clarity,  we  show  the  spectra  of  the  predicted  out¬ 
puts  (at  the  target  location)  due  to  the  SVF  (dashed-dotted), 
linear  response  (dotted),  and  the  quadratic  kernel  (dashed). 
One  can  see  that  the  linear  filter  captures  the  energy  in  the 
fundamental  band  while  the  second  order  kernel  captures 
the  energy  at  the  low  frequency  and  in  the  2nd  harmonic 
band.  The  significance  of  this  result  is  that  it  was  obtained 
based  on  training  in  the  presence  of  strong  speckle  noise 
due  to  random  scatterers  surrounding  the  target  (as  is  the 
case  in  practical  medical  imaging  applications).  Examples 


4.  CONCLUSIONS 

We  have  demonstrated  the  validity  of  a  second  order  Volterra 
system  as  a  model  for  pulse-echo  ultrasonic  imaging  data 
from  tissue  mimicking  media.  In  addition,  a  robust  method 
for  identification  of  the  SVF  coefficients  have  been  demon¬ 
strated  in  the  presence  of  speckle  noise,  an  important  result 
for  medical  ultrasound. 

With  the  recent  advances  in  VLSI  signal  processing  archi¬ 
tectures,  adaptive  imaging  algorithms  based  on  system  iden¬ 
tification  tools  are  becoming  increasingly  attractive.  The 
results  presented  in  this  paper  demonstrate  the  feasibility  of 
this  approach  to  the  emerging  class  of  nonlinear  imaging 
systems  that  will  undoubtedly  shape  the  future  of  medical 
ultrasonic  imaging. 
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Abstract-  We  describe  dual-mode  ultrasound  phased  arrays 
for  noninvaslve  image-guided  surgical  applications.  In  partic¬ 
ular,  we  address  the  problem  of  real-time  visualization  of  ther¬ 
mal  lesion  formation  In  the  target  (e.g.,  tumor)  tissue  using  the 
therapeutic  arrays.  Post  beamfroming  filter  bank  image  recon¬ 
struction  with  nonlinear  compounding  is  utilized  to  improve 
the  lesion  contrast  in  the  (typically)  very  low-contrast  ultra¬ 
sound  Images.  It  is  shown  that  the  new  image  reconstruction 
algorithm  leads  to  measurable  improvement  in  lesion  contrast 
on  the  order  of  6  -  IS  dB.  This  leads  to  significant  Improvement 
in  lesion  detectability  and  size  estimation  by  standard  segmen¬ 
tation  techniques  for  speckle  Imagery.  Experimental  results 
strongly  suggest  that  2nd  Harmonic  imaging  could  play  an  Im¬ 
portant  role  in  the  enhancement  of  real-time  lesion  visualiza¬ 
tion. 

1.  Introduction 

Image  guidance  has  long  been  recognized  as  the  “enabling 
technology”  for  noninvasive  thermal  surgery.  Highly  re¬ 
fined  energy  application  devices  have  been  developed  and 
for  years.  However,  without  reliable  imaging  techniques  for 
visualization  of  thermal  lesions,  noninvasive  thermal  surgery 
has  failed  to  find  widespread  acceptance  in  the  clinic.  Re¬ 
cently,  image  guidance  methods  based  on  well  established 
imaging  modalities  like  MRI[3],  CT[5],  or  ultrasound  [4,  8] 
have  been  proposed.  Other  imaging  modalities  are  also  be¬ 
ing  developed  and  may  become  available  in  the  foreseeable 
future[6,  7], 

An  area  unique  to  ultrasound  that  could  revolutionize 
the  field  of  image  guided  surgery  is  the  development  of  a 
new  generation  of  dual-mode  high-power  phased  array  sys¬ 
tems  capable  of  both  imaging  and  therapy  [8,  9].  These 
piezocomposite  transducers  can  produce  focal  intensity  lev¬ 
els  needed  for  ablative  and  coagulative  thermal  surgery  with 
high  precision.  Furthermore,  the  operating  bandwidth  of 
such  transducers  allows  for  imaging  the  treatment  region 
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with  adequate  image  quality  to  delineate  important  land¬ 
marks  within  and  around  the  target  volume.  With  these  ca¬ 
pabilities,  it  is  possible  to  operate  these  arrays  in  a  “self¬ 
registration”  mode  whereby  the  imaging  capabilities  of  the 
array  are  utilized  in  characterization  of  the  tissue  response 
precisely  at  the  expected  lesion  location.  This  is  due  to  the 
fact  that  the  beamforming  is  common  to  both  the  imag¬ 
ing  and  therapy  modes.  The  main  challenge  to  this  ap¬ 
proach  is  the  low-contrast  nature  of  speckle-ridden  ultra¬ 
sound  images.  This  paper  addresses  a  new  post  beamform¬ 
ing  filter  bank  image  reconstruction  algorithm  with  nonlin¬ 
ear  spatially  weighted  compounding  algorithm  to  mitigate 
this  problem.  The  design  of  this  algorithm  is  motivated  by 
the  nonlinear  nature  of  ultrasound  propagation  in  tissue  and 
the  fact  that  thermal  lesions  are  known  to  have  increased 
level  of  harmonic  generation.  Experimental  results  demon¬ 
strate  the  potential  of  the  new  algorithm  in  both  enhancing 
the  lesion  contrast  and  size/shape  analysis. 


2.  Image  Formation  Model 

Figure  1  summarizes  the  image  acquisition  and  image  for¬ 
mation  model.  A  64-element  array  optimized  for  maximum 
energy  delivery  at  1  MHz  operating  frequency  is  used  for  le¬ 
sion  formation  in  sample  tissue.  Lesions  are  formed  by  fo¬ 
cusing  the  array  at  a  point  within  the  target  and  maintaining 
high-power  output  for  time  intervals  on  the  order  of  seconds 
(1  -  5  seconds  typical).  The  power  is  interrupted  for  short 
intervals  (milliseconds)  to  acquire  image  data  by  transmit¬ 
ting  short  (ps)  pulses  from  all  64-elements  and  receiving 
on  selected  elements  using  a  matrix  switch.  Once  the  im¬ 
age  data  set  is  collected,  RF  beamforming  is  performed  to 
form  standard  echographic  images  of  the  target  region.  It 
is  well  known  that  these  images  have  very  low  contrast  due 
to  the  coherent  nature  of  image  acquisition  which  produces 
speckle.  This  renders  the  lesion  detection  very  difficult  us¬ 
ing  ultrasound. 


0-7803-7211-5/01/$17.00  ©  2001  IEEE 


2492 


*» 


Fig.  1.  Image  Acquisition  and  Image  Reconstruction  Algo¬ 
rithm 


Fig.  2.  Coordinate  system  used  in  the  synthetic  aperture 
imaging  system. 

2. 1.  Synthetic  aperture  imaging 

Images  were  obtained  using  a  full  synthetic  aperture  tech¬ 
nique  [10],  The  image  pixel  at  coordinates  ( xp,zp )  was 
therefore  computed  by  (Fig.  2): 

64  64 

I(xp,  Zp )  =  EE  Ai  •  Bj  •  Sij  [(-Rip  4-  Rjp)  /c] ,  (1) 

»=i  j=i 

where  i  is  the  transmit  element  index,  j  is  the  receive  el¬ 
ement  index,  Ai  is  the  transmit  apodization  weight  at  el¬ 
ement  i,  Bj  is  the  receive  apodization  weight  at  element 
j,  Rip  and  RjP  are  the  distances  from  the  transmit  and  re¬ 
ceive  elements,  respectively,  from  the  image  pixel,  c  is  the 
speed  of  sound  in  the  medium  being  imaged,  and  Sij(t)  is 
the  echo  acquired  when  transmitting  with  element  i  and  re¬ 
ceiving  with  element  j.  Specialized  image  reconstruction 
programs  were  written  to  evaluate  the  intensity  at  ( xp,zp ) 
in  either  (r,  sin  0)  or  Cartesian  coordinates. 


2.2.  Single-Transmit  Imaging 

The  synthetic  aperture  imaging  algorithm  described  above 
can  be  used  to  produce  the  highest  quality  conventional  im¬ 
ages  that  can  be  expected  from  a  given  array.  However,  due 
to  the  fact  that  transmit  beams  are  synthesized  by  super¬ 
position  of  single-element  transmit  patterns,  the  nonlinear 
interactions  of  the  real-time  transmit  beams  cannot  be  ac¬ 
counted  for  using  this  imaging  algorithm.  Therefore,  we 
have  modified  our  64-channel  phased  array  driver  to  allow 
for  pulsed  transmission  on  all  64  channels  simultaneously. 
This  allowed  us  to  use  the  full  power  of  the  transmit  beams 
and,  therefore,  observe  their  nonlinear  interactions  with  the 
tissue  media.  The  image  formation  process  due  to  a  single- 
transmit  focused  beam  is  a  modified  version  of  Eq.  (1)  as 
follows: 

64 

I ( xp ,  zp)  =  y  ]  Bj  •  Sj  [(i?o  4  Rjp)  /^-]  *  (2) 

i=i 

where  Sj(t)  is  the  received  waveform  at  element  j  due  to 
the  transmitted  beam  and  Ro  is  a  fixed  distance  determined 
by  the  focal  depth  of  the  transmit  beam.  All  other  quantities 
in  Eq.  (2)  are  the  same  as  their  counterparts  in  Eq.  (1). 
It  should  be  noted  that  beamforming  is  performed  in  the 
RF  domain  using  true  delays,  i.e.,  no  baseband  conversion 
is  done.  This  retains  all  the  frequency  components  in  the 
beamformed  data  for  postbeamforming  filtering  operations 
described  below. 

2.3.  Post  Beamforming  Filtering 

Since  beamformiiig  is  performed  in  the  RF  domain,  all  the 
frequency  components  in  the  received  echo  signals  are  re¬ 
tained.  This  is  important  since  the  echo  signals  can  be  ex¬ 
pected  to  contain  a  mix  of  harmonics  (and  possibly  subhar¬ 
monics)  depending  on  the  nonlinearity  of  the  tissue  medium 
being  imaged.  This  is  especially  true  for  the  single-transmit 
imaging  mode  where  the  transmit  power  is  sufficiently  high 
to  produce  harmonic  components  in  nonlinear  tissue  me¬ 
dia.  Postbeamforming  filtering  can  be  used  to  isolate  spe¬ 
cific  harmonic  components  and/or  enhance  axial  resolution 
if  the  SNR  of  the  system  is  sufficiently  high.  Algorithms  for 
postbeamforming  filterbank  image  reconstruction  for  pulse- 
echo  ultrasound  can  be  found  in  [11].  For  the  purposes  of 
this  paper,  the  filter  bank  is  formed  of  bandpass  filters  with 
center  frequencies  at  a  set  of  preselected  harmonics  (or  sub¬ 
harmonics). 

2.4.  Nonlinear  Compounding 

Ultrasound  propagation  in  tissue  media  is  inherently  nonlin¬ 
ear.  RF  echo  signals  obtained  using  modem  scanners  with 
wideband  transducers  carry  harmonic  components  that  are 
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generated  by  the  medium  nonlinearity,  i.e.,  they  are  not  part 
of  the  originally  transmitted  imaging  pulses.  Since  the  tis¬ 
sue  nonlinearity  parameter  exhibits  higher  contrast  between 
various  tissue  compared  to  tissue  reflectivity,  imaging  this 
parameter  can  lead  to  higher  contrast  images  of  soft  tissue. 
For  a  variety  of  physical  reasons,  it  is  well  known  that  ther¬ 
mal  lesions  generate  higher  harmonics  at  levels  higher  than 
normal  tissues.  By  careful  design  of  the  transmit/receive 
beamforming  and  the  post  beamforming  reconstructions  fil¬ 
ters,  high  contrast  images  of  lesions  can  be  formed  from 
ultrasound  images  formed  at  the  fundamental  and  some  of 
its  higher  harmonics.  In  general,  an  optimal  compound¬ 
ing  rule  must  be  derived  based  on  statistical  model  of  the 
imaged  region  and  the  expected  lesion  size/location.  Fur¬ 
thermore,  the  quality  of  the  beamforming  at  different  har¬ 
monics  or  scales  will  be  spatially  heterogeneous.  For  ex¬ 
ample,  since  the  array  sampling  function  (in  wavelengths) 
is  coarser  at  the  higher  harmonics,  images  formed  at  these 
frequency  will  have  high  contrast  along  the  axis  of  the  trans¬ 
mit  beam,  but  the  contrast  deteriorates  away  from  the  main 
axis.  Furthermore,  since  the  harmonics  are  orders  of  mag¬ 
nitudes  smaller  than  the  fundamental,  compounding  is  per¬ 
formed  with  log-compressed  images  after  appropriate  scal¬ 
ing.  Therefore,  with  reference  to  Figure  1,  the  compounding 
is  performed  in  two  stages: 

1.  The  outputs  of  the  filterbank  are  envelope  detected, 
scaled,  and  log  compressed. 

2.  A  spatially-weighted  sum  of  the  harmonic  compo¬ 
nents  is  performed.  The  weighting  function  is  derived 
from  the  spatial  contrast  function  (SCF)  of  the  imag¬ 
ing  array  at  the  different  harmonics. 

Assuming  that  the  envelope-detected,  normalized  and  log- 
compressed  images  obtained  from  the  filterbank  are  given 
by  Iiv,  then  the  compounded  image  is  given  by: 

N 

i(x,v)  =  '52WiSi(x’y')Ii(x’y)  <3) 

irrl 

where  Wi  is  a  weighting  function  reflecting  the  relative  en¬ 
ergy  at  the  harmonic  with  respect  to  the  fundamental  at  the 
(expected)  lesion  location  and  Si(x,y)  is  the  SCF  associ¬ 
ated  with  the  ith  harmonic.  The  SCFs  at  the  different  har¬ 
monic  components  are  computed  from  the  inverse  of  the 
beam  patterns  when  the  transmit  beam  is  focused  at  the  cen¬ 
ter  of  the  (expected)  lesion. 

3.  Results 

Images  shown  in  Figure  3  and  Figure  4  demonstrate  the 
basic  ideas  behind  the  proposed  image  compounding  algo¬ 
rithm.  The  figures  show  images  obtained  before  and  after 


lesion  formation  in  an  ex-vivo  liver  tissue  samples.  The  ex¬ 
pected  lesions  are  cigar-shaped  with  length  of  10  mm  and 
width  of  2  -  3  mm.  All  images  are  shown  at  25  dB  dynamic 
range  and  have  the  same  axial  and  lateral  extent.  In  both 
experiments,  the  lesion  was  formed  at  the  geometric  cen¬ 
ter  of  live  array  and  was  expected  between  90  and  100  mm 
in  the  axial  direction.  All  images  show  the  front  surface 
of  the  sample  at  80  mm  and  the  sponge  backing  phantom 
at  120  mm.  For  each  experiment,  two  pairs  of  images  are 
shown,  one  before  and  one  after  the  lesion  formation.  Im¬ 
ages  shown  on  the  left  hand  side  of  each  figure  show  the  tar¬ 
get  before  lesion  formation  while  images  on  the  right  hand 
side  show  the  target  after  lesion  formation.  The  images  on 
top  of  each  figure  are  formed  at  the  fundamental  (1  MHz 
50%  bandwidth)  while  those  at  the  bottom  of  each  figure 
are  formed  at  the  2nd  harmonic  (2MHz  40%  bandwidth). 
Each  pair  of  images  is  displayed  on  a  log  scale  and  normal¬ 
ized  such  that  the  0  dB  corresponds  to  the  maximum  of  the 
image  on  the  RHS,  i.e.,  the  image  of  the  target  after  lesion 
formation.  The  first  experiment  represent  a  fairly  homo- 
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Fig,  3.  Images  before  and  after  formation  of  Lesion  1 . 

geneous  target  tissue  while  the  second  has  a  pair  of  blood 
vessels  at  the  target  location  as  can  be  seen  from  the  images 
taken  before  lesion  formation.  While  both  the  fundamen¬ 
tal  and  harmonic  images  show  enhanced  lesion  contrast  in  • 

both  cases,  the  harmonic  images  show  a  net  increase  in  con¬ 
trast  by  22  dB  in  the  second  case.  On  the  other  hand,  the 
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5.  References 


Fig.  4.  Images  before  and  after  lesion  2. 


net  increase  in  contrast  at  the  fundamental  is  about  7  dB. 
In  addition,  the  spatial  definition  of  the  lesion  in  the  har¬ 
monic  images  is  superior  to  the  fundamental.  The  two  cases 
shown  also  illustrate  the  need  for  nonlinear  compounding 
with  a  spatial  weighting  function  as  implied  by  (3).  While 
the  harmonic  images  are  superior  both  in  terms  of  defini¬ 
tion  and  contrast  near  the  main  axis  of  the  transmit  beam, 
the  fundamental  images  better  represent  the  tissue  state  off- 
axis. 


4.  CONCLUSIONS 

We  have  shown  that  post  beamforming  filterbank  recon¬ 
struction  of  ultrasound  images  at  selected  frequencies  sensi¬ 
tive  to  harmonic  generation  in  nonlinear  media  produces  im¬ 
ages  appropriate  for  compounding.  The  main  features  (ther¬ 
mal  lesions)  are  sufficiently  correlated  while  the  speckle 
noise  is  largely  uncorrelated.  The  enhanced  nonlinear  tissue 
response  at  the  lesion  location  is  most  probably  due  to  for¬ 
mation  of  microbubbles  that  resonate  nonlinearly  produc¬ 
ing  2nd  and  higher  harmonic  frequencies  in  the  echo  data. 
Experimental  work  for  validation  of  this  hypothesis  is  cur¬ 
rently  underway. 
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Abstract-  A  post-beamforming  nonlinear  compounding  algo¬ 
rithm  for  ultrasonic  imaging  is  presented.  Fundamental  and 
harmonic  image  components  from  bearafromed  radio  frequency 
(RF)  data  are  extracted,  envelope  detected  and  compounded 
using  a  spatial  compounding  functions  (SCFs)  derived  from 
the  transrolt/receive  beamforming  topology  used  in  obtaining 
the  RF  data.  This  is  specially  useful  for  applications  where 
slngle-transmit  focus  (STF)  imaging  is  used,  In  this  paper,  we 
present  results  from  STF  imaging  experiments  using  a  novel 
dual-mode  phased  array  system  for  Image-guided  surgery,  In 
particular,  we  address  the  enhancement  of  the  echogenicity  of 
thermal  lesions  formed  in  ex  vivo  tissue.  It  is  shown  that  new 
nonlinear  image  compounding  algorithm  produces  25  -  30  dB 
enhancement  In  lesion  echogenicity  without  loss  in  spatial  res¬ 
olution.  This  Is  to  be  compared  with  a  typical  enhancement 
of5  dB  achieved  by  standard  echographlc  imaging  and  20  dB 
achieved  by  second  harmonic  imaging  alone.  In  addition,  im¬ 
ages  resulting  from  the  new  algorithm  are  virtually  Bee  of 
beamforming  artifacts  that  can  severely  degrade  the  perfor¬ 
mance  of  2nd  harmonic  imaging. 

1,  INTRODUCTION 

Image  guidance  has  long  been  recognized  as  the  “enabling 
technology"  for  noninvasive  thermal  surgery.  Highly  re¬ 
fined  energy  application  devices  have  been  developed  and 
for  years.  However,  without  reliable  imaging  techniques  for 
visualization  of  thermal  lesions,  noninvasive  thermal  surgery 
has  failed  to  find  widespread  acceptance  in  the  clinic.  Re¬ 
cently,  image  guidance  methods  based  on  well  established 
imaging  modalities  like  MRI[3],  CT[5),  or  ultrasound  [4,  8] 
have  been  proposed.  Other  imaging  modalities  are  also  be¬ 
ing  developed  and  may  become  available  in  the  foreseeable 
futureI6, 7]. 

An  area  unique  to  ultrasound  that  could  revolutionize 
the  field  of  image-guided  surgery  is  the  development  of  a 
new  generation  of  dual-mode  high-power  phased  array  sys- 
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terns  capable  of  both  imaging  and  therapy  [8,  9],  These 
piezocomposite  transducers  can  produce  focal  intensity  lev¬ 
els  needed  for  ablative  and  coagulative  thermal  surgery  with 
high  precision.  Furthermore,  the  operating  bandwidth  of 
such  transducers  allows  for  imaging  the  treatment  region 
with  adequate  image  quality  to  delineate  important  land¬ 
marks  within  and  around  the  target  volume.  With  these  ca¬ 
pabilities,  it  is  possible  to  operate  these  arrays  in  a  “self- 
registration”  mode  whereby  the  imaging  capabilities  of  the 
array  are  utilized  in  characterization  of  the  tissue  response 
precisely  at  the  expected  lesion  location.  This  is  due  to  the 
fact  that  the  beamforming  is  common  to  both  the  imag¬ 
ing  and  therapy  modes.  The  main  challenge  to  this  ap¬ 
proach  is  the  low-contrast  nature  of  speckle-ridden  ultra¬ 
sound  images.  This  paper  addresses  a  new  post  beamform¬ 
ing  filter  bank  image  reconstruction  algorithm  with  nonlin¬ 
ear  spatially  weighted  compounding  algorithm  to  mitigate 
this  problem.  The  design  of  this  algorithm  is  motivated  by 
the  nonlinear  nature  of  ultrasound  propagation  in  tissue  and 
the  fact  that  thermal  lesions  are  known  to  have  increased 
level  of  harmonic  generation.  Experimental  results  demon¬ 
strate  the  potential  of  the  new  algorithm  in  both  enhancing 
the  lesion  contrast  and  size/shape  analysis. 


2.  IMAGE  FORMATION  MODEL 

Figure  l  summarizes  the  image  acquisition  and  image  for¬ 
mation  model.  A  64-eiement  array  optimized  for  maximum 
energy  delivery  at  1  MHz  operating  frequency  is  used  for  le¬ 
sion  formation  in  sample  tissue.  Lesions  are  formed  by  fo¬ 
cusing  the  array  at  a  point  within  the  target  and  maintaining 
high-power  output  for  time  intervals  on  the  order  of  seconds 
(1-5  seconds  typical).  The  power  is  interrupted  for  short 
intervals  (milliseconds)  to  acquire  image  data  by  transmit¬ 
ting  short  (ps)  pulses  from  all  64-eiements  and  receiving  on 
selected  elements  using  a  matrix  switch.  Once  the  image 
data  set  is  collected,  RF  beamforming  is  performed  to  form 
standard  echographic  images  of  the  target  region. 
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Fig.  1.  Image  Acquisition  and  Image  Reconstruction  Algorithm 


Fig.  2.  Coordinate  system  used  in  the  synthetic  aperture 
imaging  system. 

2.1.  Synthetic  aperture  imaging 

Images  are  obtained  using  a  full  synthetic  aperture  tech¬ 
nique  [10],  The  image  pixel  at  coordinates  (a?p,  zp)  is  there¬ 
fore  computed  by  (Fig.  2): 

61  64 

zp)  =  At '  '  st,i  +  Rjp)  /c]  >  (1) 

t=t  j= i 

where  *  is  the  transmit  element  index,  j  is  the  receive  el¬ 
ement  index,  Aj  is  the  transmit  apodization  weight  at  ele¬ 
ment  t,  Bj  is  the  receive  apodization  weight  at  element  j. 
Rip  and  Rjp  are  the  distances  from  the  transmit  and  receive 
elements,  respectively,  from  the  image  pixel,  c  is  the  speed 
of  sound  in  the  medium  being  imaged,  and  s,j(t)  is  the 
echo  acquired  when  transmitting  with  element  i  and  receiv¬ 
ing  with  element  j. 

2.2.  Single-Transmit  Focus  (STF)  Imaging 

The  synthetic  aperture  imaging  algorithm  described  above 
can  be  used  to  produce  the  highest  quality  conventional  im¬ 
ages  that  can  he  expected  from  a  given  array.  However,  due 
to  the  fact  that  transmit  beams  are  synthesized  by  super¬ 
position  of  single-element  transmit  patterns,  the  nonlinear 
interactions  of  the  real-time  transmit  beams  cannot  be  ac¬ 
counted  for  using  this  imaging  algorithm.  Therefore,  we 


have  modified  our  64-channel  phased  array  driver  to  allow 
for  pulsed  transmission  on  all  64  channels  simultaneously. 
This  allowed  us  to  use  the  full  power  of  the  transmit  beams 
and,  therefore,  observe  their  nonlinear  interactions  with  the 
tissue  media.  The  image  formation  process  due  to  a  STF 
beam  is  a  modified  version  of  Equation  (1)  as  follows: 

61 

J(xp,  Zp)  =  ^  Bj  •  sj  [{Ro  +  Rjp)  / c] ,  (2) 

J=i 

where  Sj  (t)  is  the  received  waveform  at  element  j  due  to  the 
transmitted  beam  and  i?o  is  a  fixed  distance  determined  by 
the  focal  depth  of  the  transmit  beam.  All  other  quantities  in 
Equation  (2)  are  the  same  as  their  counterparts  in  Equation 
(I).  It  should  be  noted  that  beamforming  is  performed  in  the 
RF  domain  using  true  delays,  i.e.,  no  baseband  conversion 
is  done.  This  retains  all  the  frequency  components  in  the 
beamformed  data  for  postbeamforming  filtering  operations 
described  below. 

23.  Post  Beamforming  Filtering 
Since  beamforming  is  performed  in  the  RF  domain,  all  the 
frequency  components  in  the  received  echo  signals  are  re¬ 
tained.  This  is  important  since  the  echo  signals  can  be  ex¬ 
pected  to  contain  a  mix  of  harmonics  (and  possibly  subhar¬ 
monics)  depending  on  the  nonlinearity  of  the  tissue  medium 
being  imaged.  This  is  especially  true  for  the  STF  imaging 
mode  where  the  transmit  power  is  sufficiently  high  to  pro¬ 
duce  harmonic  components  in  nonlinear  tissue  media.  Post¬ 
beamforming  filtering  can  be  used  to  isolate  specific  har¬ 
monic  components  and/or  enhance  axial  resolution  if  the 
SNR  of  the  system  is  sufficiently  high.  Algorithms  for  post- 
beamforming  interbank  image  reconstruction  for  pulse-echo 
ultrasound  can  be  found  in  [11],  For  the  purposes  of  this 
paper,  the  filter  bank  is  formed  of  bandpass  filters  with  cen¬ 
ter  frequencies  at  a  set  of  preselected  harmonics  (or  subhar- 
monics). 

2.4.  Nonlinear  Compounding 

Ultrasound  propagation  in  tissue  media  is  inherently  nonlin¬ 
ear.  RF  echo  signals  obtained  using  modem  scanners  with 
wideband  transducers  carry  harmonic  components  that  are 
generated  by  the  medium  nonlinearity,  i.e.,  they  are  not  part 
of  the  originally  transmitted  imaging  pulses.  Since  the  tis¬ 
sue  nonlinearity  parameter  exhibits  higher  contrast  between 
various  tissue  compared  to  tissue  reflectivity,  imaging  this 
parameter  can  lead  to  higher  contrast  images  of  soft  tissue. 
For  a  variety  of  physical  reasons,  it  is  well  known  that  ther¬ 
mal  lesions  generate  higher  harmonics  at  levels  higher  than 
normal  tissues.  By  careful  design  of  the  transmit/receive 
beamforming  and  the  post  beamforming  reconstructions  fil¬ 
ters,  high  contrast  images  of  lesions  can  he  formed  from 
ultrasound  images  formed  at  the  fundamental  and  some  of 
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its  higher  harmonics.  In  genera],  an  optimal  compound¬ 
ing  rule  must  be  derived  based  on  statistical  model  of  the 
imaged  region  and  the  expected  lesion  size/location.  Fur¬ 
thermore,  the  quality  of  the  beamforming  at  different  har¬ 
monics  or  scales  will  be  spatially  heterogeneous.  For  exam¬ 
ple,  since  the  array  sampling  function  (in  wavelengths)  is 
coarser  at  the  higher  harmonics,  images  formed  at  these  fre¬ 
quency  will  have  high  contrast  along  the  axis  of  the  trans¬ 
mit  beam,  but  the  contrast  deteriorates  quickly  away  from 
the  main  axis.  Furthermore,  since  the  harmonics  are  orders 
of  magnitudes  smaller  than  the  fundamental,  compounding 
is  performed  with  log-compressed  images  after  appropriate 
scaling.  Therefore,  with  reference  to  Figure  1,  the  com¬ 
pounding  is  performed  in  two  stages: 

1.  The  outputs  of  the  filterbank  are  envelope  detected, 
scaled,  and  log  compressed.  The  scaling  of  the  indi¬ 
vidual  harmonics  can  be  derived  from,  for  example, 
the  specificity  or  contrast-to-tissue  ratio  (CTR): 

CTR  =  1°'  0‘,"(mI)  (3) 

where  ||yc,.  H2  and  ||yT(  ||2  are  the  norms  of  the  tth 
harmonic  components  from  the  contrast  and  normal 
tissue  regions,  respectively.  These  regions  are  easily 
identified  under  various  imaging  conditions.  For  in¬ 
stance,  for  the  application  described  in  this  paper,  the 
contrast  region  is  the  expected  location  of  the  ther¬ 
mal  lesion  (often  visible  on  the  standard  echographic 
image). 

2.  A  spatially-weighted  sum  of  the  harmonic  compo¬ 
nents  is  performed.  The  weighting  function  is  derived 
from  the  spatial  contrast  function  (SCF)  of  the  imag¬ 
ing  array  at  the  different  harmonics. 

Assuming  that  the  envelope-detected,  normalized  and  log- 
compressed  images  obtained  from  the  filterbank  are  given 
by  h,  I2, . . . ,  In,  then  the  compounded  image  is  given  by: 

N 

I{x>  y)  Wi$i(x, y)Ii{x,  y)  (4) 
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where  Wt  is  a  weighting  function  reflecting  the  relative  en¬ 
ergy  at  the  harmonic  with  respect  to  the  fundamental  at  the 
(expected)  lesion  location  and  $i(x,y )  is  the  SCF  associ¬ 
ated  with  the  ith  harmonic.  This  is  can  be  obtained  based 
on  simulations  of  pulsed  or  CW  calculations  of  the  follow¬ 
ing  ratio: 


Si(x,  y) 


Jml  b<-<  (*>  V'~  x')Bt(*,  Voi  x'W 
Jsl  Bn  (*.  VI  x')Bt(x:  j/o ;  x')dx’  ' 


where  the  subscripts  si  and  ml  refer  to  sidelobe  and  main- 
lobe,  respectively,  and  Bri  (x,  y,  x1)  is  the  tth  harmonic  (dy¬ 
namic)  receive  focus  at  (x,  y)  and  Bt(x,  yo;x')  is  the  trans¬ 
mit  focus  pattern  (at  the  fundamental).  Note  that  the  SCF 


is  a  measure  of  the  quality  of  the  transmit/receive  pattern  at 
the  tth  harmonic.  The  higher  the  SCF  the  more  confidence 
we  have  in  the  received  data  at  the  tth  harmonic. 

3.  RESULTS 

Figure  3  gives  a  typical  result  of  a  visualization  experiment 
using  the  dual-mode  array  and  the  image  compounding  al¬ 
gorithm  described  above.  The  figure  shows  images  obtained 
before  and  after  lesion  formation  in  an  ex- vivo  liver  tissue 
samples.  The  expected  lesions  are  cigar-shaped  with  length 
of  10  mm  and  width  of  2  -  3  mm.  All  images  are  shown 
at  40  dB  dynamic  range  and  have  the  same  axial  and  lateral 
extent  The  lesion  was  formed  at  the  geometric  center  of  the 
array  and  was  expected  between  90  and  100  mm  in  the  axial 
direction.  All  images  show  the  front  surface  of  the  sample  at 
80  mm  and  the  sponge  backing  phantom  at  120  mm.  Three 
pairs  of  images  are  shown,  one  before  (left)  and  one  after 
the  lesion  formation.  The  images  on  top  of  each  figure  are 
formed  at  the  fundamental  (1  MHz)  while  those  at  the  center 
of  each  figure  are  formed  at  the  2nd  harmonic  (2  MHz).  The 
pair  at  the  bottom  shows  the  compound  image  using  Equa¬ 
tion  4  and  the  SCFs  shown  in  Figure  4.  Each  pair  of  images 
is  displayed  on  a  log  scale  and  normalized  such  that  the  0 
dB  corresponds  to  the  maximum  of  the  image  on  the  RHS, 
i.e„  the  image  of  the  target  after  lesion  formation.  While 
both  the  fundamental  and  harmonic  images  show  enhanced 
lesion  contrast,  the  harmonic  images  show  a  net  increase  in 
contrast  by  22  dB  along  the  direction  of  the  transmit  beam. 
On  the  other  hand,  the  net  increase  in  contrast  at  the  funda¬ 
mental  is  about  7  dB.  In  addition,  the  spatial  definition  of 
the  lesion  in  the  harmonic  images  is  superior  to  the  funda¬ 
mental.  However,  while  the  harmonic  images  are  superior 
both  in  terms  of  definition  and  contrast  near  the  main  axis  of 
the  transmit  beam,  the  fundamental  images  better  represent 
the  tissue  state  off-axis.  The  compound  images  show  a  net 
increase  in  the  lesion  contrast  on  the  order  of  35  dB  both 
axially  and  laterally.  This  result  can  be  understood  in  light 
of  a  close  examination  of  the  SCRs  shown  in  Figure  4,  It  is 
also  interesting  to  note,  however,  that  the  compound  image 
still  shows  major  off-axis  objects,  e.g.,  blood  vessel  at  axial 
distance  of  1 10  mm  and  lateral  distance  of  5  mm.  That  is, 
the  compounding  algorithm  removes  the  beamforming  side- 
lobe  artifacts,  but  not  major  scatterers  in  the  vicinity  of  the 
transmit  beam.  This  is  one  of  the  main  objectives  of  STF 
imaging  described  in  this  paper. 


4.  CONCLUSIONS 

We  have  shown  that  post  beamforming  filterbank  recon¬ 
struction  of  ultrasound  images  at  selected  frequencies  sen¬ 
sitive  to  harmonic  generation  in  nonlinear  media  produces 
images  appropriate  for  compounding.  The  results  of  apply¬ 
ing  the  nonlinear  compounding  algorithm  also  demonstrate 
a  significant  increase  in  lesion  contrast  over  enhancement 
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Fig.  3.  Images  before  (left)  and  after  lesion  formation.  All  images 
are  shown  at  40  dB  dynamic  range. 


achieved  at  either  the  fundamental  or  second  harmonic.  It  is 
important  to  state  that  this  contrast  enhancement  is  achieved 
without  loss  in  spatial  resolution.  This  due  to  the  fact  that 
the  bandpass  filters  used  in  separating  the  fundamental  and 
harmonic  components  have  wider  bandwidth  than  the  data 
bandwidth  at  these  harmonics.  The  results  shown  in  this  pa¬ 
per  suggest  that  dual-mode  arrays  with  high  contrast  imag¬ 
ing  capability  suitable  for  real-time  thermal  lesion  visual¬ 
ization  are  feasible.  This  leads  to  a  most  powerful  paradigm 
for  image  guided  surgery  where  the  therapeutic  and  image- 
guidance  coordinate  systems  are  inherently  registered. 
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Abstract  -  This  paper  provides  a  first  experimental  verifica¬ 
tion  of  what  is  arguably  the  most  significant  benefit  of  dual-mode 
ultrasound  array  (DMUA)  systems  for  noninvasive  surgery.  Specif¬ 
ically,  the  use  of  image-based  feedback  for  refocusing  the  therapeu¬ 
tic  beam  in  the  presence  of  strongly  scattering  objects.  This  capa¬ 
bility  may  be  critical  in  the  use  of  DMUA  systems  for  noninvasive 
targeting  of  liver  tumors  and  noninvasive  cardiac  ablation.  In  both 
cases,  the  target  is  partially  obstructed  by  the  rib  cage,  which  lim¬ 
its  the  access  and  distorts  the  geometrically-focused  high-intensity 
focused  ultrasound  (HIFU)  therapeutic  beam.  The  imaging  capa¬ 
bilities  of  the  applicator  are  used  in  the  optimization  of  the  tar¬ 
geting  process.  The  optimization  procedure  is  based  on  the  use 
of  single-transmit  focus  (STF)  imaging,  in  which  a  single  trans¬ 
mit  imaging  beam  employing  delays  derived  from  phasing  for  the 
therapeutic  beam  is  used.  Full  dynamic  receive  focusing  is  used 
to  form  images  of  the  target  region.  Taking  advantage  of  the  in¬ 
herent  registration  between  the  imaging  and  therapeutic  coordinate 
systems,  STF  images  can  be  used  in  refocusing  the  DMUA  at  the 
target(s)  while  minimizing  the  power  deposition  at  the  rib(s)  within 
the  HIFU  beam.  We  present  experimental  results  of  a  DMUA  re¬ 
focusing  algorithm  that  selects  the  control  points  from  the  targetfs) 
and  the  ribs  visible  in  the  STF  image.  The  new  phases  and  ampli¬ 
tudes  (for  HIFU  beam)  or  focusing  delays  (for  adaptive  focusing) 
are  obtained  by  solving  a  constrained  optimization  problem  that 
minimizes  the  direct  incidence  of  acoustic  power  at  the  ribs  as  it 
optimizes  the  power  delivered  to  the  target.  A  64-element  1  MHz 
DMUA  prototype  was  used  to  refocus  on  a  fine  needle  thermocou¬ 
ple  within  a  tissue-mimicking  phantom  near  the  geometric  focus. 
Four  Plexiglas  1.2  cm  diameter  ribs  were  placed  in  front  of  the 
target  (as  seen  from  the  DMUA  prototype).  Thermocouples  were 
positioned  on  the  ribs  on  the  DMUA  side  in  the  imaging/therapy 
plane.  The  temperature  at  the  target  and  the  ribs  were  recorded  be¬ 
fore,  during  and  after  3  s  HIFU  exposure  for  the  following  cases: 
1)  geometric  focusing,  2)  refocusing  without  considering  the  pres¬ 
ence  of  the  ribs,  and  3)  refocusing  with  accounting  for  the  ribs.  The 
input  electrical  power  to  the  system  was  kept  practically  the  same 
and  the  ratios  of  temperature  at  the  target  and  the  ribs  per  watt  of 
input  power  were  recorded.  The  results  of  these  experiments  show 
that  the  power  delivered  at  the  target  can  be  increased  by  several 
fold  (due  to  refocusing)  without  a  corresponding  increase  at  the 
rib  locations.  This  conclusion  is  supported  by  direct  thermocouple 
measurement  at  the  target  and  rib  locations.  In  addition,  STF  im¬ 
ages  before  and  after  refocusing  provide  similar  feedback  by  con¬ 
sistently  showing  increased  echogenicity  of  the  target  region  while 
the  echogenicity  of  the  ribs  is  not  increased  (often  reduced) 


I.  Introduction 

Ultrasound  phased  array  systems  are  especially  suitable 
for  the  generation  of  HIFU  beams  in  the  presence  of  tis¬ 
sue  inhomogeneities  and/or  strong  obstacles  such  as  bones 
and  skull.  Recent  results  from  Fink  and  co-workers  [1]  and 
Hynynen  and  co-workers  [2]  have  shown,  for  example,  the 
feasibility  of  generating  HIFU  beams  through  the  skull  with 
proper  phase  feedback.  Both  groups  proposed  refocusing 
based  on  time  reversal  (phase  conjugation)  by  means  of  a 
receiving  hydrophone  at  the  focus  location.  They  have  also 
proposed  a  computational  approach  based  on  an  acoustic 
model  derived  from  3D  CT  or  MRI  image  of  the  target  vol¬ 
ume.  We  have  also  investigated  the  optimal  synthesis  of 
HIFU  beams  through  the  rib  cage  using  a  two-step  approach 
to  minimize  the  incident  power  on  the  ribs  while  maximiz¬ 
ing  the  power  delivered  to  the  target  (e.g.  tumor)  [3], 

The  introduction  of  piezocomposite  transducer  technol¬ 
ogy  in  the  fabrication  of  therapeutic  ultrasound  phased  ar¬ 
rays  has  created  a  new  opportunity  to  investigate  a  new 
generation  of  dual-mode  ultrasound  array  (DMUA)  systems 
for  imaging  and  therapy  [4].  We  have  recently  demon¬ 
strated  both  the  therapeutic  and  imaging  capabilities  of  a  64- 
element  1  MHz  array  prototype.  HIFU-induced  lesions  in  ex 
vivo  tissue  demonstrated  a  high  degree  of  exposure  control 
within  a  therapeutic  operating  field  ( ThxOF )  extending  by 
±1.5  cm  around  the  geometric  center  of  the  array.  With  a 
37%  fractional  bandwidth,  it  was  shown  that  the  array  has 
an  imaging  field  of  view  ( IxFOV )  extending  well  beyond 
the  ThxOF  [5].  Experimental  results  in  the  visualization 
of  HIFU-induced  lesions  in  ex  vivo  tissue  have  shown  that 
DMUA  images  allow  for  mapping  of  the  lesion  boundaries 
based  on  changes  in  echogenicity  [6].  Have  shown  a  far 
greater  level  of  consistency  than  those  obtained  with  stan¬ 
dard  diagnostic  equipment,  which  was  somewhat  surprising. 
However,  this  may  be  attributed  to  two  important  reasons:  1) 
imaging  with  the  same  array  allows  the  interrogation  of  the 
exact  treatment  site  due  to  the  fact  that  the  imaging  and  ther¬ 
apeutic  beams  are  aligned,  and  2)  imaging  with  the  same 
frequency  of  the  therapeutic  beam  is  more  sensitive  to  any 
resonant  bubble  activity  that  may  have  been  generated  by 
rectified  diffusion. 

In  addition  to  the  above  advantage  of  DMUA  with  regards 
to  lesion  visualization,  the  imaging  capabilities  of  these  ar¬ 
rays  offer  another  unique  advantage  with  regards  to  guid¬ 
ance  of  HIFU  beams.  In  particular,  DMUA  images  can  pro- 
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vide  feedback  for  refocusing  the  HIFU  beam  taking  into  ac¬ 
count  any  obstacles  or  critical  structures  that  may  be  inter¬ 
acting  with  the  beam.  One  application  of  this  capability  is 
refocusing  in  the  presence  of  the  ribs  when  targeting  liver  tu¬ 
mors  or  in  noninvasive  cardiac  ablation.  We  have  formulated 
this  as  an  optimal  synthesis  of  phased  array  beam  patterns 
with  control  points  at  the  target  and  the  critical  structures  (in 
this  case  the  ribs).  The  optimization  procedure  is  based  on 
the  use  of  single-transmit  focus  (STF)  imaging,  in  which  a 
single  transmit  imaging  beam  employing  the  same  delays  as 
the  intended  therapeutic  focus  is  used.  Full  dynamic  receive 
focusing  is  used  to  form  images  of  the  target  region.  The  re¬ 
sulting  image  is  sensitive  to  the  STF  beam  intensity  and  the 
strength  of  the  scattering  structure  along  the  beam  path.  It 
is  expected  that  the  ribs  will  produce  strong  echoes  as  long 
as  they  are  within  the  main  beam  or  significant  sidelobes. 
Therefore,  they  are  easily  identified  in  the  STF  images  al¬ 
lowing  for  the  estimation  of  the  DMUA  element  directivity 
at  these  locations.  The  target  point(s),  on  the  other  hand, 
can  be  located  on  the  STF  image  at  desired  refocus  loca- 
tion(s),  i.e.  no  need  for  strong  scatterer  at  target  location. 
With  DMUA  element  directivity  at  both  target  and  obsta¬ 
cles  defined  from  beamforming  delays  used  to  form  the  STF 
image,  a  constrained  optimization  problem  maximizing  the 
power  deposition  at  the  target  while  minimizing  the  power 
deposition  at  the  ribs  can  be  formulated.  In  what  follows,  we 
give  a  description  for  DMUAs  image  formation,  constrained 
optimization  refocusing  algorithm,  and  experimental  verifi¬ 
cation  of  this  algorithm. 

II.  THEORY  AND  METHODS 

A.  Synthetic  Aperture  Imaging 

Images  were  obtained  using  a  full  synthetic  aperture  tech¬ 
nique  [7].  This  was  also  described  in  our  previous  publica¬ 
tions  [5].  This  gives  the  highest  quality  image  possible  with 
the  DMUA  and  allows  for  most  accurate  visualization  of  the 
treatment  region. 

B.  Single-Transtnit  Imaging 

We  have  proposed  the  use  of  this  imaging  mode  for  vi¬ 
sualization  of  discrete  lesion  formation  in  previous  reports 
[5].  However,  one  very  important  application  of  this  mode 
is  in  surveying  the  target  region  for  any  critical  structures  in¬ 
teracting  with  the  therapeutic  beam.  By  using  beamforming 
delays  for  the  imaging  transmit  beams  corresponding  to  the 
phases  of  the  therapeutic  beams,  a  dynamic  receive  image  is 
a  result  of  the  relative  strength  of  the  beam  and  any  scatterers 
within  the  main  beam  and  significant  sidelobes.  For  exam¬ 
ple,  strongly  scattering  obstacles  such  as  the  ribs  within  the 
beam  are  always  visible  on  the  STF  image  of  the  DMUA. 
Since  the  beamforming  (delays)  for  any  point  on  the  image 
are  known,  the  locations  of  these  obstacles  can  be  used  as 
control  points  in  synthesizing  refocused  therapeutic  beams 


that  minimize  the  energy  at  the  ribs  while  maximizing  the 
energy  at  the  target. 

C.  Refocusing  in  the  Presence  of  Obstacles 

1)  Refocusing  Algorithm:  Assume  we  have  an  In¬ 
clement  array  with  arbitrary,  but  known,  geometry  radi¬ 
ating  into  a  homogeneous  half  space.  The  objective  is 
to  maximize  the  array  intensity  gain  at  a  target  point  fr 
while  minimize  the  incident  power  at  a  set  of  critical 
points,  tq  (i),  i  =  1,2,  Me-  This  is  an  optimization 
problem  that  can  be  solved  using  Lagrange  multipliers  or  a 
regularized  minimum-norm  least  squares  solution  as  follows 
[8]: 

•  Let  the  row  vector  hj  = 

[^i(Ur),  h^fr),  ■  ■  ■  ,  ^At(rr)]  be  the  array  direc¬ 
tivity  vector  at  i~t  (note  that  h-kirr)  is  the  directivity 
of  the  kth  element  at  fp)- 

•  Form  the  matrix  He  from  the  row  vectors  h  = 
[hi(fc(i)),h2(fc(i)),---  ,hN(rc(i))\. 

•  Form  the  weighting  matrix,  Wq  =  [HcH£f  +  ^yl] — 
where  I  is  the  identify  matrix  and  7  is  an  appropriately 
chosen  regularization  parameter. 

•  The  optimal  complex  array  driving  vector  is  given  by 
the  u  =  Wc’hljr  (hj-Wc’hJf  )-1. 

Both  the  target  and  the  critical  points  can  be  derived  from 
images  formed  by  the  DMUA  in  synthetic  aperture  or  in 
single-transmit  focus  modes.  This  is  probably  the  most  pow¬ 
erful  feature  of  the  use  of  the  DMUA  in  image-guided  appli¬ 
cation. 

2)  Refocusing  Procedure:  Based  on  the  above  algorithm, 
the  refocusing  procedure  is  as  follows: 

•  SA  and/or  STF  images  are  formed  to  survey  the  target 
and  any  strongly  scattering  objects  in  the  therapeutic 
beam  path  (or  significant  sidelobes). 

•  The  elements  of  hj-  and  He  are  determined  from  the 
images  and  the  optimal  refocused  vector  is  evaluated. 

•  The  optimal  steering  vector  is  used  to  obtain  an  STF 
image  for  assessment  of  the  refocusing  result.  If  the  re¬ 
focusing  is  successful,  then  the  target  region  should  ap¬ 
pear  brighter  while  the  critical  structures  should  appear 
less  bright  compared  to  the  image  before  refocusing. 

•  If  STF  image  feedback  indicates  successful  refocusing, 
the  optimal  steering  vector  is  used  therapeutically  and 
the  direct  effects  of  heating  are  measured  at  the  target 
and  critical  structures,  e.g.  temperature  measurements. 

III.  Experiment  Setup 

Figure  1  shows  the  image  acquisition  model.  A  64- 
element  concave  linear  phased  array  (100  mm  radius  of  cur¬ 
vature)  is  used  for  the  image  acquisition.  For  synthetic  aper¬ 
ture  images,  a  single-cycle  1  MHz  pulse  is  transmitted  by 
each  individual  element  and  received  by  each  of  the  ele¬ 
ments  to  generate  the  RF  data.  Once  the  RF  data  set  is 
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Figure  1 :  image  experiment  setup 
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Figure  2:  schematic  diagram  of  the  phantom 

collected,  RF  beamforming  is  performed  to  form  standard 
echographic  images  of  the  target  region.  For  STF  images, 
a  single  transmit  beam  for  all  elements  is  used  and  the  echo 
data  is  received  on  individual  elements,  i.e.  64  transmissions 
are  needed  to  obtain  a  complete  STF  data  set. 

A.  Tissue-Mimicking  phantom  with  obstacles 

In  order  to  demonstrate  the  feasibility  of  optimization 
of  power  deposition  at  the  target(s)  while  minimizing  the 
power  deposition  at  strongly  scattering  obstacles,  we  de¬ 
signed  the  experiment  setup  shown  in  Figure  2.  Four  8- 
mm  diameter  plexiglas  ribs  with  spacing  of  20  mm  are  used 
to  simulate  the  ribcage.  They  are  positioned  about  50  mm 
away  from  the  dual-mode  array.  A  tissue-mimicking  phan¬ 
tom  [9][10]  (75mmx60  mmx37mm)  is  positioned  near  the 
geometric  focus  of  the  DMUA  (behind  the  ribs).  A  Needle 
thermocouple  (D=0.2  mm)  is  placed  near  the  geometric  fo¬ 
cus  (not  necessarily  at  the  geometric  focus).  Another  two 
wire  thermocouples  are  fixed  on  the  two  middle  ribs  in  the 
ribs  cage.  The  three  thermocouples  are  measured  every  0.1 
second  by  the  HP  34970A  data  acquisition  unit. 

IV.  Results 

A.  Imaging:  Suryery,  Refocusing,  and  Assessment 

Figure  3  shows  a  grayscale  image  (50  dB)  of  the  phan¬ 
tom.  The  image  was  obtained  using  the  full  synthetic  aper¬ 
ture  technique  described  in  [5].  The  result  is  consistent  with 
the  schematic  with  good  spatial  resolution  as  well  as  con¬ 
trast.  In  addition  to  the  target  and  the  ribs,  the  front  and  back 
of  the  tissue-mimicking  phantom  can  be  easily  discerned. 


Figure  3:  Grayscale  image  of  the  tissue-mimicking  phan¬ 
tom  with  ribs  shown  in  Figure  2  acquired  using  the  DMA  in 
synthetic  aperture  mode  (50  dB).  The  white  marks  to  left  are 
1  cm  apart. 


Even  though  the  ribs  may  be  outside  the  50  dB  IxFOV  of  the 
DMUA,  they  are  sufficiently  strong  to  standout  giving  strong 
specular  reflections.  Furthermore,  the  ribs  are  all  spatially 
registered  correctly  without  distortion.  The  image  shown  in 
Figure  3  was  used  to  identify  the  coordinates  of  the  visible 
ribs  as  well  as  a  target  point  for  refocusing  within  the  tissue- 
mimicking  phantom.  The  40  dB  grayscale  images  shown  in 
Figure  4  show  the  results  of  refocusing  before  (left)  and  af¬ 
ter  (right)  taking  the  rib  location  into  consideration.  Both  of 
these  images  where  obtained  with  a  transmit  beam  focused 
at  the  defined  target  location  at  -4  mm  lateral  and  103  mm 
axial.  One  can  clearly  see  the  increased  visibility  of  the  tar¬ 
get  compared  to  the  visibility  of  the  ribs.  As  shown  is  Figure 
5,  after  refocusing,  the  echo  from  the  target  area  is  enhanced 
by  10  dB  while  at  the  rib  area,  the  echo  is  reduced  by  7  dB. 
This  is  due  to  the  reduction  of  incident  power  at  the  ribs  as 
predicted  by  the  weighting  algorithm  described  in  II-C.l. 

B.  Direct  Temperature  Measurement 

Figure  6  shows  the  temperature  changes  at  the  target  and 
two  rib  locations  during  a  3-second  exposure  using  thera¬ 
peutic  beams  corresponding  to  the  imaging  beams  used  for 
generating  the  STF  images  in  Figure  4,  i.e.  before  and  after 
refocusing.  The  temperature  is  normalized  by  the  DC  power 
delivered  to  the  DMUA.  Clearly,  the  incident  power  to  the 
target  is  increased  by  a  factor  of  5  to  6,  while  the  power  at 
the  rib  locations  is  increased  by  a  factor  of  two  or  less. 

V.  Conclusions 

This  paper  presented  the  first  experimental  verification 
of  the  feasibility  of  using  image-based  feedback  for  refo¬ 
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(a)  (b) 


Figure  4:  STF  images  (40  dB)  using  refocusing  before  (a) 
and  after  (b)  taking  the  critical  points  and  the  ribs  into  con¬ 
sideration. 


(a)  (b) 

Figure  5:  The  magnitude  change  in  dB  in  (a)  lateral  direc¬ 
tion  passing  by  the  target;(b)  axial  direction  passing  by  the 
rib. 

cusing  for  targeting  tissue  structures  while  avoiding  critical 
obstacles  such  as  the  ribs.  This  may  be  important  if  HIFU 
arrays  are  to  be  used  for  targeting,  for  example,  liver  tu¬ 
mors  that  may  be  partially  obstructed  by  rib  cages.  This 
result  may  not  be  understood  as  a  form  of  aberration  cor¬ 
rection  since  the  refocusing  algorithm  assumes  knowledge 
of  the  array  directivity  at  the  refocusing  target  and  the  crit¬ 
ical  points.  Nonetheless,  since  these  same  assumptions  are 
used  in  beamforming  the  image  in  the  first  place,  the  algo¬ 
rithm  is  quite  robust  against  distortion  due  to  inaccuracies 
in  the  speed  of  sound  and  other  tissue  properties.  The  only 
requirement  is  that  at  the  SA  or  the  STF  image  used  in  refo- 


Figure  6:  Temperature  changes  at  the  ribs  and  target  lo¬ 
cation  without  and  with  the  refocusing  algorithm  from:  (a) 
target;  (b)  rib  1;  (c)  rib  2. 


cusing  provides  a  recognizable  map  of  the  treatment  region. 
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Noninvasive  Localized  Ultrasonic  Measurement  of  Tissue  properties 
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Abstract —  We  present  in  vitro  and  in  vivo  results  for  valida¬ 
tion  of  a  new  localized  noninvasive  ultrasonic  measurement  of  both 
tissue  absorption  and  perfusion.  The  method  employs  sub-second 
low-intensity  focused  ultrasonic  beams  for  generating  brief  temper¬ 
ature  rise  on  the  order  of  1°  C.  The  RF  data  from  an  imaging  trans¬ 
ducer  is  processed  to  produce  noninvasive  temperature  estimate  in 
the  localized  heated  volume  using  speckle  tracking  or  frequency- 
domain  algorithms  previously  published  [1]  [2],  Absorption  can  be 
obtained  from  the  initial  heating  rate  while  perfusion  can  be  esti¬ 
mated  from  the  initial  decay.  Due  to  the  small  size  of  the  heated 
spot,  the  noninvasive  temperature  estimation  can  produce  results 
that  are  virtually  free  of  thermal  lensing  effects.  Furthermore,  since 
the  measurements  are  based  on  heating  and  decay  rates,  only  the 
parameters  of  the  transient  bioheat  equation  (density  and  heat  ca¬ 
pacity)  are  needed  for  the  estimate.  Both  ex  -vivo  and  in  vivo  results 
have  shown  two-fold  to  three-fold  increase  in  tissue  absorption.  For 
the  ex  vivo  results,  the  change  in  absorption  was  estimated  using  di¬ 
rect  fine  wire  thermocouple  measurements  at  the  treatment  site  in 
addition  to  the  noninvasive  temperature  estimation.  The  decay  rate 
of  the  in  vivo  estimated  temperature  was  observed  to  increase  by 
two-fold  indicating  increased  perfusion  in  the  tumor  surrounding 
the  small  lesion.  While  the  opposite  effect  can  be  expected  in  volu¬ 
metric  lesion  formation,  this  is  very  likely  in  this  single-shot  lesion 
formation  experiment.  In  any  event,  the  in  vivo  results  show  very 
clearly  the  feasibility  of  estimating  perfusion  based  on  decay  rate. 

I.  Introduction 

High  intensity  focused  ultrasound  (HIFU)  is  a 
promising  modality  for  image-guided  noninvasive 
surgery.  MRI  is  currently  the  leading  modality  for 
guidance  of  HIFU  due  to  its  excellent  quality  and  tem¬ 
perature  sensitivity.  Ultrasound  guidance  of  HIFU  is 
limited  by  the  non-quantitative  nature  of  conventional 
ultrasonic  imaging.  Several  groups  have  recently  inves¬ 
tigated  the  use  of  elastography  and  viscoelastic  remote 
palpation  techniques  for  characterization  of  the  state  of 
ablated  tissue  (both  HIFU  and  RF  ablation)  [3].  Ini¬ 
tial  reports  on  remote  palpation  are  promising  but  true 
localization  of  the  shear  waves  to  allow  a  simple  esti- 
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mation  of  the  shear  elastic  modulus  is  yet  to  be  demon¬ 
strated  in  tissue. 

Two  quantities  believed  to  reflect  the  state  of  the  tis¬ 
sue  after  HIFU  lesion  formation  are  tissue  absorption 
and  blood  perfusion.  The  former  has  been  shown  to 
increase  upon  lesion  formation  while  the  latter  is  hy¬ 
pothesized  to  exhibit  dynamics  depending  on  the  tem¬ 
perature  rise,  position  within  the  heated  volume,  and 
degree  of  damage  to  the  vasculature.  The  proposed 
new  tissue  property  measurement  is  based  on  nonin¬ 
vasive  tissue  temperature  estimation  due  to  a  short  ex¬ 
posure  (on  the  order  of  10-1  s)  from  a  focused  beam  to 
create  temperature  rise  on  the  order  of  1°C.  This  expo¬ 
sure  is  therapeutically  insignificant  and  may  fall  well 
within  the  range  of  allowable  thermal  indices  of  stan¬ 
dard  imaging  equipment.  A  variety  of  geometries  for 
the  heating  source  and  the  imaging  probe  can  be  con¬ 
sidered.  However,  in  this  paper,  we  present  results  from 
a  setup  where  the  imaging  plane  is  orthogonal  to  the 
axis  of  the  heating  transducer.  The  temperature  estima¬ 
tion  method  used  in  this  paper  is  based  on  the  thermal 
dependence  of  the  ultrasound  echo  that  account  for  two 
different  physical  phenomena:  local  changes  of  sound 
speed  due  to  temperature  variation  and  thermal  expan¬ 
sion  of  the  propagating  medium.  The  former  generates 
an  apparent  shift  in  scattered  location,  and  the  latter 
leads  to  a  physical  shift.  Along  an  A-line,  however,  the 
two  effects  lead  to  echo  time-shift  that  can  be  estimated 
and  have  been  shown  to  be  related  local  change  in  tem¬ 
perature  in  the  propagating  medium.  This  temperature 
estimation  approach  is  based  on  a  method  originally 
described  in  [1]  and  a  later  development  given  in  [2]. 

II.  Localized  Tissue  Property  Estimation 

The  2D  temperature  estimation  follows  the  algo¬ 
rithm  given  in  [2]  and  is  summarized  here  to  illustrate 
its  use  for  the  purposes  of  the  parameter  estimation  al¬ 
gorithm  (see  [2]  for  full  detail).  The  cumulative  echo 
time-shift  5t(z,  x,  I) )  at  each  location  and  time  T,;,  are 
estimated  by  a  correlation-based  speckle  tracking  algo- 
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rithm.  Axial  differentiation  on  these  estimates  (along 
the  ^  direction)  is  performed  in  conjunction  with  2D 
spatial  filtering  along  both  axial  ( z )  and  the  lateral  (x) 
direction.  The  symbol  T)  corresponds  to  the  wall  clock 
time,  at  which  the  i-th  frame  was  acquired  and  it  should 
be  distinguished  from  echo  time-delays  t.  The  algo¬ 
rithm  can  be  described  by  the  following  steps: 

Stepl  i  —  0.  Acquire  2-D  RF-data  r(z,  x,  To),  prior 
to  any  heating  (baseline  temperature  9q). 

Step2  i  —  i  +  1.  Acquire  2-D  RF-data  r(z,  x,T:) 
every  AT  0.1. 

If  Tinit  <  Ti  <  Tend,  Heating  Beam  is  ON. 

Step3  Estimate  the  incremental  time-shift  map 


Technos  MPX 

1 

system 

fitincr(z,  Tj)  —  X,Ti)  i[z,Xi,T%—\) 

at  time  T)  using  the  current  and  previous 
frames  r(z,  x,  Ti)  and  r(z,  x,  T,_i). 

Step4  Compute  the  cumulative  time-shift  map 

i 

^iincr{,tt,X,Ti)  —  ^  fiijncr(.Z,  X,  Tfc). 

k= 1 

Step5  Differentiate  the  cumulative  time-shift  map 
S6(z,  x,  Tj)  at  time  Ti  (including  lowpass  fil¬ 
tering  in  spatial  and  temporal  directions). 

Step6  If  Ti  <  Tmox,  Go  to  Step2. 

In  this  algorithm,  Tinit,Tend,Tmax,  and  AT  denote 
start  heating  time,  stop  heating  time,  end  of  acquisition 
time,  and  frame  period,  respectively.  Parameter  estima¬ 
tion  is  performed  using  the  transient  temperature  pro¬ 
file  at  the  maximum  spatial  peak  in  the  2D  temperature 
rise.  The  absorption  is  proportional  to  the  initial  heat¬ 
ing  rate  while  perfusion  is  related  to  the  initial  decay 
rate  upon  turning  off  the  heating  beam.  These  param¬ 
eters  are  obtained  assuming  the  transient  bioheat  equa¬ 
tion  (TBHTE)  model,  which  is  widely  accepted  in  the 
hyperthermia  community  [4]. 

IIT.  Experiment  Setup 

All  imaging  was  performed  using  a  modified  Tech- 
nos  MPX  system  from  ESAOTE  S.p.A.,  Genoa,  Italy. 
The  Technos  is  equipped  with  an  RF  grabber  for  cap¬ 
turing  high  quality  beamformed  RF  data,  which  al¬ 
lowed  the  capture  of  up  to  60  seconds  of  full  frame  data 
with  a  specified  frame  rate.  A  linear  array  probe  (LA 
522,  1-cycle  transmit  pulse  centered  at  8.5  MHz)  was 


Figure  1:  Experiment  setup  used  for  HIFU-induced 
lesion  formation  using  a  single  element  transducer  and 
observation  using  an  LA  522  linear  probe  monitoring 
along  the  axial  direction  of  the  lesion 


heating  beam 

Figure  2:  Experiment  time  schedule 

used  in  acquiring  the  RF  data  during  the  heating  experi¬ 
ments.  A  4MHz  transducer  with  40  mm  focus  was  used 
for  in  vivo  for  low-temperature  short-duration  heating 
experiments  (Focus  Surgery,  Inc.).  A  1.5  MHz  spher¬ 
ical  shell  transducer  with  63.5  mm  radius  of  curvature 
was  used  for  ex  vivo  and  phantom  experiments. 

Each  data  set  was  25  seconds  of  frame  data  (at  4 
frame/s)  from  the  Technos  acquired  as  described  in 
Section  II.  The  first  3  seconds  consisted  of  frames  be¬ 
fore  heating  with  the  focused  transducer,  i.e.  Tmax  — 
25  s  and  Tinu  —  3  s.  Then  the  HIFU  beam  was  turned 
on  for  0.2  second  or  0.5  second  as  shown  in  Figure  2. 
Note  that  the  actual  measurement  time  is  on  the  order 
of  1  s  and  not  the  total  25  s  used  to  collect  the  data. 

A.  Experiment  Materials 

Experiments  were  performed  on  tissue  mimicking 
phantom  material  as  described  in  [5][3],  degassed  ex 
vivo  porcine  liver,  and  one  nu/nu  immunodeficient 
mouse  in  vivo.  The  mouse  had  an  MA148  human  ovar¬ 
ian  carcinoma  grown  in  the  hindleg.  The  elastography 
phantoms  were  fabricated  from  gelatin,  graphite,  alco¬ 
hol,  water  and  glutaraldehyde  as  described  in  the  cita- 
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Figure  3:  (a)2-D  estimated  temperature  change  im¬ 

mediately  after  turning  off  the  focused  transducer; 
Normalized  temperature  changes  of  tissue-mimicking 
phantom  under  different  exposure  conditions  from 
thermocouple  (b)  and  temperature  estimation  algo¬ 
rithm  (c).  (d)  Comparison  of  initial  heating  rate  from 
noninvasive  temperature  estimation  with  thermocouple 
measurements.  The  maximum  temperature  change  is 
4.2  °C  and  the  mimimum  temperature  change  is  0.75 
°C. 


tions. 


IV.  Results 

A.  Tissue-mimicking  Phantom 

Figure  3  shows  the  measured  temperature  changes 
collected  from  a  needle  thermocouple  (D=0.2  nun) 
near  the  transducer  focus  location  and  the  noninva¬ 
sive  temperature  estimation  results  at  the  focus  location 
from  the  ultrasound  echoes.  The  estimated  tempera¬ 
ture  profiles  have  similar  transients  as  the  thermocou¬ 
ple  measurements.  The  results  clearly  show  excellent 
agreement  between  the  estimated  and  measured  results 
in  terms  of  the  initial  heating  rate  and  decay  rate.  Fig¬ 
ure  3  (d)  shows  that  1°C  change  for  less  than  0.5  s  can 
be  tracked  very  well  using  the  noninvasive  temperature 
estimation  technique. 

To  further  quantify  the  agreement  between  the  esti¬ 
mated  and  measured  temperatures,  we  estimate  the  ini¬ 
tial  decay  rate  immediately  after  powering  off  the  heat¬ 
ing  transducer  (within  1  s  of  power  off).  Table  I  shows 


TABLE  I 

Initial  Decay  Rate 


Exposure 

Noninv.  Est. 

Thermocouple 

1W  0.2s 

-0.169 

-0.116 

1W  0.5s 

-0.111 

-0.115 

4W  0.2s 

-0.099 

-0.105 

4W  0.5s 

-0.084 

-0.103 

the  initial  decay  rate  for  the  various  curves  shown  in 
Figure  3  (b)  and  (c).  The  results  from  the  two  sets 
of  curves  are  consistent,  except  for  the  lowest  expo¬ 
sure  level  resulting  in  0.75°  temperature  rise.  This  is 
also  consistent  with  Figure  3(d),  which  allows  the  same 
conclusion. 

B.  In-vitro  Experiment  Result 

Figure  4  shows  the  measured  temperature  changes 
collected  from  a  needle  thermocouple  near  the  trans¬ 
ducer  focus  location  and  the  temperature  estimation  re¬ 
sults  at  the  focus  location  from  the  ultrasound  echoes 
gotten  from  a  degassed  ex  vivo  porcine  liver.  Initial 
heating  rates  from  the  estimated  and  direct  measure¬ 
ments  are  in  good  agreement  before  and  after  lesion 
formation.  An  increase  in  the  initial  heating  rate  at  the 
lesion  location  by  a  factor  of  2.5  is  observed  here. 

C.  In-vivo  Experiment  Result 

Figure  5  shows  the  estimated  transient  temperature 
profiles  at  the  focus  location  before  and  after  HIFU  le¬ 
sion  formation  in  MAI 48  human  ovarian  carcinoma 
in  the  hindleg  of  the  nu/nu  immunodeficient  mouse. 
Comparison  between  the  profiles  before  and  after  le¬ 
sion  formation  indicates  that  a  factor  of  2.5  increase  in 
the  initial  heating  rate  has  been  achieved  in  this  case 
too.  In  addition,  it  appears  that  the  decay  rate  after  le¬ 
sion  formation  has  increased.  While  this  is  the  result 
of  only  one  experiment,  it  may  be  due  to  an  increase 
in  perfusion  in  the  vasculature  surrounding  the  small 
single-shot  lesion  that  was  formed  in  this  case.  This  is 
the  subject  of  an  ongoing  research  in  our  laboratory. 

V.  Conclusions 

The  experimental  results  shown  in  this  paper  demon¬ 
strate  the  feasibility  of  reliable  estimation  of  the  ini¬ 
tial  heating  rate  at  a  localized  heating  spot  for  short- 
duration  small  temperature  changes  on  the  order  of  1° 
in  tissue  media,  including  in  vivo.  Furthermore,  the  ini¬ 
tial  decay  rates  can  be  reliably  measured  in  phantoms 
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(a)  (b) 


Figure  4:  Normalized  temperature  changes  of  pig  liver 
under  different  conditions  from:  (a)  thermocouple  be¬ 
fore  lesion  formation;  (b)  thermocouple  after  lesion 
formation;  (c)  estimated  results  before  lesion  forma¬ 
tion;  (d)  estimated  results  after  lesion  formation.  The 
maximum  temperature  change  is  1.2  °C  and  the  mini¬ 
mum  temperature  change  is  0.2°C. 

as  well  as  tissue  media  for  the  same  type  of  low  ex¬ 
posure  heating.  Ongoing  research  at  our  laboratory  in¬ 
vestigates  the  feasibility  of  measuring  local  tissue  per¬ 
fusion  from  the  initial  decay  rate  of  the  localized  tissue 
heating.  The  2D  temperature  estimation  results  indi¬ 
cate  that  the  spatial  resolution  of  such  measurements  is 
determined  by  the  heating  source  and  the  imaging  sys¬ 
tem.  For  the  configuration  used  for  this  paper,  2  mm 
axial  resolution  is  possible  (see  Figure  3(a)). 
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Real-Time  Monitoring  of  the  Transients  of  HIFU-Induced  Lesions 


Hui  Yao  and  Emad  S.  Ebbini 
Department  of  Electrical  and  Computer  Engineering 
University  of  Minnesota,  Minneapolis,  MN  55455 


Abstract  -  RF  data  from  standard  B-mode  and  pulse 
inversion  (PI)  imaging  of  FIIFU  lesion  formation  of  freshly 
excised  tissue  was  collected  before,  during,  and  after  lesion 
formation  experiment's  in  ex  vivo  tissue.  Exposures  at  inten¬ 
sity  levels  of  1100  W/cm2  to  2500W/cm2  for  durations  of  2, 
3,  and  5  seconds  in  a  single  shot  were  used.  Also  continuous 
raster  scan  of  longer  duration  (10  -  20  seconds)  to  form  slice 
or  volumetric  lesions  were  monitored.  Monitoring  was  done 
with  a  diagnostic  scanner  and  RF  data  was  acquired  at  1 
frame/second  for  60  seconds  starting  5  seconds  before  each 
shot.  Lesion  maps  from  grayscale  B-mode  and  PI  images 
were  obtained  using  level-set  methods  for  each  frame  and 
compared  with  the  actual  lesion  found  by  inspection.  Le¬ 
sion  maps  from  PI  imaging  were  consistently  smaller  in  size 
and  more  in  line  with  the  actual  lesion  size.  Transient  anal¬ 
ysis  of  harmonic  content  of  lesion  echoes  show  sustained 
harmonic  activity  for  1 0  - 1 5  seconds  after  the  therapy  pulse 
is  turned  off  (in  ex  vivo  liver  tissue).  A  gradual  drop  in  this 
activity  follows  with  steady  state  reached  within  50  -  60  sec¬ 
onds.  It  was  also  shown  that  the  use  of  short  microsecond 
pulses  from  the  therapy  transducer  to  expose  the  lesion  lo¬ 
cation  during  real-time  imaging  significantly  increased  the 
scattering  from  lesion  location. 

I.  Introduction 

Ultrasound  is  currently  a  leading  imaging  modalities  for 
guidance  of  numerous  interventional  procedures.  This  is  due 
to  advantages  in  cost,  portability,  ease  of  integration,  and 
real-time  implementation.  However,  its  compromised  con¬ 
trast  resolution  and  poorly  understood  HIFU  lesion  dynam¬ 
ics  present  major  obstacles  to  the  goal  of  mapping  these  le¬ 
sions  with  conventional  ultrasound.  Over  the  last  few  years, 
several  research  groups  suggested  the  use  of  more  quantita¬ 
tive  techniques  such  as  elastography  and  remote  palpation 
for  mapping  HIFU  lesions.  The  successful  implementation 
of  these  methods  relies  on  better  understanding  of  the  tran¬ 
sients  of  the  echoes  from  the  HIFU  lesion  location. 

Motivated  by  P.  P.  Lele  in  [1],  we  hypothesized  that  this 
change  in  echogenicity  is  due  to  stable  microbubbles  that 
can  occur  even  at  low  insonation  levels.  Lele  found  that 
subharmonic  emission  due  to  microbubbles  showed  a  mono¬ 
tonic  increase  with  intensity  from  150  mW/cm2  to  1500 
W/cm2  without  a  distinct  threshold  for  emission  (measure¬ 
ment  done  in  vitro  and  in  vivo  at  2.7  and  1.8  MHz).  The 


consistency  of  the  increase  in  echogenicity  at  the  lesion  may 
be  explained  by  the  fact  that  the  microbubbles  may  already 
be  resonant  at  the  imaging  frequency  (same  as  the  therapeu¬ 
tic  HIFU  beam). 

II,  Theory  and  Mathematical  Models 

A.  Level  set  theory 

Reference  [2]  gives  a  multi  phase  level  set  framework  for 
image  segmentation  using  the  Mumford  and  Shah  model  for 
piecewise  constant  and  piecewise  smooth  optimal  approx¬ 
imation.  This  level  set  method  is  especially  suitable  for 
HlFU-induced  lesion  size  estimation  due  to  their  irregular 
shape. 

Given  an  observed  image  uq,  to  find  a  decomposition  Q; 
of  H  and  an  optimal  piecewise  smooth  approximation  u  of 
Uo  such  that  u  varies  smoothly  inside  Q,  and  discontinu¬ 
ous^  across  the  boundaries  of  Q,  is  called  the  segmentation 
problem  in  computer  vision  [3j. 

Given  uo,  find  (u,  C)  to  minimize 

irf{EMS(u,C)}  (1) 

where, 

Ems(u,C)  =  f  ( u-uo)2dxdy+p,  [  \Vu\2dxdy+v\C\ 
Jrt  Jn\c 

(2) 

where  p,  u  >  0  are  fixed  parameters. 

A  reduced  case  of  the  above  model  is  called  minimal  par¬ 
tition  problem : 

Ems[ ti,  C)  =  f  ( u  -  uo)2 dxdy  +  v\C\  (3) 

J  n( 

In  the  piecewise-constant  case,  two  level  set  functions 
produce  four  phases  (classes)  as  shown  in  Figure  1 . 

III.  Experiment  Setup 

Figure  2  shows  a  basic  experiment  setup  for  the  formation 
of  HIFU-induced  lesions  in  freshly  excised  and  degassed 
porcine  livers  samples.  Experiments  were  performed  with 
a  modified  Technos  MPX  system  (ESAOTE  SPA,  Genoa, 
Italy)  and  a  therapy  transducer  (Piezo  Technologies,  Etalon, 
Indianapolis,  IN).  The  Technos  system  was  modified  to  al¬ 
low  imaging  in  pulse  inversion  mode  in  addition  to  normal 
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Figure  1 :  Illustration  of  two  level  set  functions:  FiSure  3:  0ne  programmed  movement  of  the  phantom 

holder  to  form  a  volumetric  lesion. 


ninr 

Inal  rtp*« 


Figure  2:  Experiment  setup  used  for  HIFU-induced  lesion 
formation  using  a  single  element  transducer  and  observation 
using  CA  421  imaging  probe  monitoring  along  the  axial  di¬ 
rection  of  the  lesion. 


B-mode  imaging.  In  addition,  a  hardware  module  for  cap¬ 
turing  high  quality  beamformed  RF  data  allows  us  to  cap¬ 
ture  and  upload  up  to  60  seconds  of  full  frame  data  with  a 
specified  frame  rate.  Frame  and  line  trigger  signals  from  the 
Technos  system  are  also  available  for  synchronizing  exter¬ 
nal  events  with  image  acquisition.  For  this  paper,  image  data 
were  acquired  from  a  CA  421  convex  probe  in  PI  mode  with 
a  2-cycle  transmit  pulse  centered  at  1 .85  MHz.  The  therapy 
transducer  is  a  1 .5  MHz  single-element  spherical-shell  trans¬ 
ducer  with  /#  =  1  (and  focal  distance  of  63.5  mm.)  A  HP 
function  generator  33120  A  (Hewlett  Packard)  is  connected 
to  the  therapy  transducer  through  an  ENI  Model  A 150  RF 
power  Amplifier  as  the  H1FU  source. 

A.  Single  shot  lesion  formation 

A  variety  of  exposures  conditions,  before,  during,  and  af¬ 
ter  lesion  formation  were  tested.  RF  data  from  standard  B- 
mode  and  PI  imaging  of  HIFU  lesion  formation  of  freshly 
excised  tissue  was  obtained.  Exposures  at  intensity  levels 
of  II 00  W/cm2  to  2500  W/cm2  for  durations  of  2,  3,  and  5 
seconds  in  a  single  shot  were  used. 


B.  Volumetric  lesion  formation 

Volumetric  lesions  are  formed  by  moving  the  tissue  holder 
at  a  constant  speed  while  the  therapeutic  beam  is  on  with 
an  intensity  suitable  for  lesion  formation.  Figure  3  gives 
one  example  of  a  raster  scan  used  in  forming  a  10  x  8  x  10 
mm3  lesion.  The  translation  speed  has  to  be  carefully  set  to 
produce  the  desired  lesion  shape  (without  over  exposure)  for 
a  given  therapeutic  intensity. 

C.  Imaging  with  the  therapeutic  beam 

Our  previous  work  on  imaging  with  dual-mode  arrays 
(DMAs)  ex  vivo  has  shown  that  the  echogenicity  of  HIFU- 
induced  lesions  is  consistently  higher,  even  minutes  after  the 
lesion  has  formed  [4],  However,  the  same  is  not  true  when 
monitoring  lesions  with  standard  diagnostic  imaging  trans¬ 
ducers.  In  order  to  gain  a  better  understanding  of  this  is¬ 
sue,  we  have  designed  a  synchronized  imaging  mode  with 
the  modified  Technos  MPX  scanner  that  can  be  described  as 
follows: 

1 .  The  setup  of  Figure  2  is  used. 

2.  A  line  trigger  signal  is  obtained  from  the  Technos 
scanner,  which  is  running  at  an  extremely  low  me¬ 
chanical  index  (MI)  in  PI  mode. 

3.  On  each  line  trigger,  the  HP  function  generator 
sends  a  2  cycle  1.5  MHz  sinusoidal  pulse  with  24 
ps  delay  through  the  therapeutic  transducer.  The 
pulse  from  the  therapeutic  transducer  has  an  ex¬ 
tremely  low  MI  as  well  (0.038  as  measured  by  a 
calibrated  hydrophone).  The  delay  is  computed 
so  that  scattering  from  the  lesion  location  registers 
correctly  on  the  Technos  image. 

4.  The  RF  data  is  uploaded  and  PI  images  are  formed. 
Since  the  imaging  pulses  are  so  small,  almost  per¬ 
fect  cancellation  of  the  backscattering  is  achieved 
throughout  the  image.  However,  since  the  pulses 
from  the  therapeutic  transducer  have  the  same  po¬ 
larity  for  both  the  positive  and  negative  PI  pulses, 
they  add  up  upon  formation  of  the  PI  image.  Thus, 
the  image  contains  the  angular  scattering  compo¬ 
nents  from  the  therapeutic  beam. 

5.  Thirty  frames  (1  fps)  are  acquired  before  lesion 
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Figure  4:  Second  harmonic  and  ultra-harmonic  energy  as 
a  function  of  time  starting  immediately  after  the  therapeutic 
beam  is  turned  off. 

formation,  1 5  without  and  1 5  with  the  pulse  from 
the  therapeutic  beam. 

6.  Fifty  ftames  (1  fps)  are  acquired  immediately  after 
lesion  formation,  28  without  and  22  with  the  pulse 
from  the  therapeutic  beam. 

This  experiment  is  intended  as  an  approximation  of  our 
DMA  system.  Here  we  are  using  the  Technos  to  image 
scattering  from  the  therapeutic  beam  at  very  low  intensity 
and  microsecond  pulse  duration.  In  the  DMA  approach,  the 
same  thing  is  done  to  obtain  the  backscatter  from  the  thera¬ 
peutic  beam. 

IV.  RESULTS 

A.  Transient  harmonic  content  of  lesion 

Figure  4  illustrate  a  typical  transient  analysis  result  of  har¬ 
monic  content  of  lesion  echoes  after  the  therapy  pulse  was 
-turned  off.  The  data  is  typical  of  an  over  exposure  single¬ 
shot  lesion  (2500  W/cm2  for  5  seconds.)  Strong  harmonic 
activity  appears  to  persist  for  10  -  15  s  at  the  second  har¬ 
monic  (2  /o)  and  the  ultra  harmonic  (1 .5  /0)  before  decaying 
to  steady  state  values  as  can  be  seen. 

B.  Lesion  mapping 

Lesion  maps  from  grayscale  B-mode  and  PI  images  were 
obtained  using  level-set  methods  for  each  frame  and  com¬ 
pared  with  the  actual  lesion  found  by  inspection.  Figure  5  il¬ 
lustrates  volumetric  lesion  maps  from  grayscale  B-mode  and 
PI  images  obtained  using  level-set  method.  The  volumetric 
lesion  was  formed  by  moving  the  phantom  holder  at  2  mm/s 
while  maintaining  the  intensity  levels  at  2500  W/cm2.  The 
lesion  size  in  B-mode  images  is  172.81  mm2,  in  PI  images 
is  1 1 1.19  mm2  while  the  real  lesion  size  is  83.08  mm2.  Le¬ 
sion  maps  from  PI  imaging  were  consistently  smaller  in  size 
and  more  in  line  with  the  actual  lesion  size.  Furthermore, 
the  transient  behavior  of  the  lesion  maps  was  easier  to  ob¬ 
serve  when  the  contours  were  evaluated  using  the  PI  images. 
Please  note  that  the  lesion  maps  given  here  are  meant  to  de¬ 
fine  the  boundaries  of  the  actual  lesion.  We  simply  want  to 
identify  the  tissue  region,  as  seen  on  an  ultrasound  image, 
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Figure  5:  Cross  section  of  a  volumetric  lesion  in  ex  vivo 
porcine  liver:  (a)  Picture;  (b)  Lesion  map  from  grayscale  B- 
mode  obtained  by  the  level-set  method;  (c)  Lesion  map  from 
grayscale  PI  mode  obtained  by  the  level-set  method. 

that  showed  a  change  in  echogenicity  due  to  the  formation 
of  a  lesion  by  a  HIFU  beam.  Quantitative  assessment  of  ir¬ 
reversible  change  of  tissue  due  to  lesion  formation  is  not  an 
objective  of  this  paper. 

C.  Imaging  with  therapeutic  beams 

Figure  6  demonstrates  the  value  of  synchronized  imag¬ 
ing  with  microsecond  low  MI  pulses  from  the  therapeutic 
transducer  as  described  in  Section  III-C.  The  spectrograms 
(SPGs)  show  the  frequency  content  of  the  received  RP  signal 
along  an  A-line  through  the  lesion  location.  The  top  SPGs 
are  computed  using  data  in  normal  B-mode  imaging.  Fig¬ 
ure  6(a),  (b),  and  (c)  are  SPGs  from  frames  acquired  before, 
1  s  after,  and  28  s  after  lesion  formation.  One  can  see  a  sig¬ 
nificant  increase  in  harmonic  activity  after  lesion  formation, 
including  sub-harmonic  activity  1  s  after  lesion  formation. 
The  sub-harmonic  component  is  not  visible  28  second  af¬ 
ter  lesion  formation.  Furthermore,  the  extent  of  the  region 
where  the  harmonic  activity  occurs  is  also  diminished  28  s 
after  lesion  formation.  This  is  consistent  with  the  results 
from  Section  IV-A.  Figure  6  (d),  (e),  and  (f)  are  SPGs  ob¬ 
tained  before,  29  s  after,  and  5 1  s  after  lesion  formation, 
all  synchronized  with  the  2  cycle  l  .5  MHz  pulses  from  the 
therapeutic  transducer  as  described  in  Section  III-C.  One 
can  see  that  the  presence  of  the  pulse  from  the  therapeutic 
transducer  before  lesion  formation  only  shows  as  a  1 .5  MHz 
component  at  the  lesion  location,  but  no  other  nonlinearity 
is  observed.  On  the  other  hand,  the  SPGs  29  and  51  s  af¬ 
ter  lesion  formation  are  very  similar  in  terms  of  the  signif¬ 
icant  harmonic  activity,  including  sub-harmonic  generation. 
A  possible  explanation  of  this  result  is  that  even  very  low 
MI  pulses  at  the  frequency  of  the  therapeutic  beam  are  capa¬ 
ble  of  generating  significant  harmonic  activity  at  the  lesion 
location.  It  is  possible  that  microbubbles  generated  by  the 
HIFU  beam  still  exist  well  after  the  transients  have  died  and 
that  even  microsecond  pulses  at  their  resonance  frequency 
are  capable  of  producing  significant  increase  in  both  linear 
and  nonlinear  scattering  from  the  lesion  location. 

Figure  7  shows  the  energy  from  the  lesion  location  be¬ 
fore  (bottom)  and  after  (top)  lesion  formation.  The  energy  is 
computed  from  received  RF  data  as  described  in  Section  HI- 
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V.  CONCLUSIONS 


Figure  6:  Spectrograms  (SPGs)  of  the  RF  signal  from  the 
center  line  passing  through  the  lesion.  Top:  w/o  (a)  before, 
(b)  1  s,  and  (c)  28  s  after.  Bottom:  w/  (d)  before,  (e)  29  s, 
and  (f)  51s  after. 


Figure  7:  Total  energy  from  the  lesion  location  (a)  after 
lesion  formation;  (b)  before  lesion  formation.  Before  the 
lesion  was  formed,  the  therapy  beam  caused  about  25  dB 
energy  increase.  The  break  in  the  curve  corresponds  to  the 
introduction  of  the  2  cycle  1 .5  MHz  pulse  from  the  thera¬ 
peutic  transducer  as  described  in  Section  III-C. 


C  first  without  and  later  with  the  synchronized  2  cycle  1.5 
MHz  pulses  from  the  therapeutic  transducer.  Without  the 
insonation  using  the  therapeutic  transducer,  the  energy  ap¬ 
pears  to  decay  after  lesion  formation  in  a  similar  manner  to 
what  we  described  in  Section  IV-A.  However,  once  the  in¬ 
sonation  using  the  therapeutic  transducer  begins,  a  sustained 
increase  in  the  energy  is  achieved,  33  dB  above  the  base  line 
in  this  case.  Examining  the  behavior  of  the  energy  before 
lesion  formation,  the  increase  due  to  the  therapeutic  beam 
is  only  25  dB.  This  indicates  that  the  scattering  from  the 
lesion  location  is  significantly  higher  after  lesion  formation. 
This  is  more  in  line  with  our  previously  reported  results  with 
the  DMA  and  illustrates  the  inherent  advantage  of  imaging 
HIFU  induced  lesions  at  the  same  frequency  as  the  therapeu¬ 
tic  beams. 


The  results  presented  in  this  paper  show  that  the  increase 
in  echogenicity  at  the  HIFU  lesion  location  can  be  transient 
in  nature  when  imaged  using  standard  diagnostic  scanners. 
In  addition,  the  size  of  the  hyperechoic  region  in  B-mode 
images  can  vary  greatly  and  may  not  be  well  correlated  with 
the  actual  lesion  size  in  some  cases.  Quantitative  measure¬ 
ments  of  tissue  property  that  may  indicate  irreversible  tissue 
damage  immediately  upon  lesion  formation  is  not  very  likely 
based  on  these  standard  images. 

On  the  other  hand,  nonlinear  (e.g.  pulse  inversion)  imag¬ 
ing  techniques  consistently  provide  tighter  maps  of  lesion 
location  on  the  ultrasound  that  is  typically  better  correlated 
with  tissue  damage.  This  mode  of  imaging,  in  addition  to 
probing  the  lesion  site  with  the  same  therapeutic  beam  (at 
low  Ml)  may  very  well  provide  the  tools  for  addressing 
this  important  problem  in  ultrasound  monitoring  of  HIFU- 
induced  lesions.  The  results  shown  herein  further  support 
our  working  hypothesis  that  imaging  HIFU-induced  lesions 
with  the  therapeutic  beam  using  dual-mode  arrays  may  be  a 
key  component  in  the  future  of  image-guided  HIFU  systems 
for  noninvasive  surgery. 
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